Function to fit selection models.

selmodel(x, ...)

# S3 method for rma.uni
selmodel(x, type, alternative="greater", prec,
         delta, steps, verbose=FALSE, digits, control, ...)

Arguments

x

an object of class "rma.uni".

type

character string to specify the type of selection model. Possible options are "beta", "halfnorm", "negexp", "logistic", "power", "negexppow", or "stepfun". Can be abbreviated. See ‘Details’.

alternative

character string to specify the sidedness of the hypothesis when testing the observed effect sizes or outcomes. Possible options are "greater" (the default), "less", or "two.sided". Can be abbreviated.

prec

optional character string to specify the measure of precision (only relevant for selection models that can incorporate this into the selection function). Possible options are "sei", "vi", "ninv", or "sqrtninv". See ‘Details’.

delta

optional numeric vector (of the same length as the number of selection model parameters) to fix the corresponding \(\delta\) value(s). A specific \(\delta\) value can be fixed by setting the corresponding element of this argument to the desired value. A specific \(\delta\) value will be estimated if the corresponding element is set equal to NA. See ‘Details’.

steps

numeric vector of one or more values between 0 and 1 that can or must be specified for certain selection functions. See ‘Details’.

verbose

logical to specify whether output should be generated on the progress of the model fitting (the default is FALSE). Can also be an integer. Values > 1 generate more verbose output. See ‘Note’.

digits

integer to specify the number of decimal places to which the printed results should be rounded (if unspecified, the default is to take the value from the object).

control

optional list of control values for the estimation algorithm. If unspecified, default values are defined inside the function. See ‘Note’.

...

other arguments.

Details

Selection models are a general class of models that attempt to model the process by which the studies included in a meta-analysis may have been influenced by some form of publication bias. If a particular selection model is an adequate approximation for the underlying selection process, then the model provides estimates of the parameters of interest (e.g., the average true effect or outcome and the amount of heterogeneity in the true effects or outcomes) that are ‘corrected’ for this selection process (i.e., they are estimates of the parameters in the population of studies before any selection has taken place). The present function fits a variety of such selection models. To do so, one should pass an object fitted with the rma.uni function to the first argument. The model that will then be fitted is of the same form as the original model combined with the specific selection model chosen (see below for possible options). For example, if the original model was a random-effects model, then a random-effects selection model will be fitted. Similarly, if the original model included moderators, then they will also be included in the selection model. Model fitting is done via maximum likelihood (ML) estimation over the fixed- and random-effects parameters (e.g., \(\mu\) and \(\tau^2\) in a random-effects model) and the selection model parameters.

Argument type determines the specific type of selection model that should be fitted. All selection models that can be fitted are based on the idea that selection may haven taken place based on the p-values of the studies. In particular, let \(y_i\) and \(v_i\) denote the observed effect size or outcome and the corresponding sampling variance of the \(i\)th study. Then \(z_i = y_i / \sqrt{v_i}\) is the (Wald-type) test statistic for testing the null hypothesis \(\mbox{H}_0{:}\; \theta_i = 0\) and \(p_i = 1 - \Phi(z_i)\) (if alternative="greater"), \(p_i = \Phi(z_i)\) (if alternative="less"), or \(p_i = 2(1 - \Phi(|z_i|)\) (if alternative="two.sided") the corresponding (one- or two-sided) p-value, where \(\Phi()\) denotes the cumulative distribution function of a standard normal distribution. Finally, let \(w(p_i)\) denote some function that specifies the relative likelihood of selection given the p-value of a study.

If \(w(p_i) > w(p_{i'})\) when \(p_i < p_{i'}\) (i.e., \(w(p_i)\) is larger for smaller p-values), then alternative="greater" implies selection in favor of increasingly significant positive effects, alternative="less" implies selection in favor of increasingly significant negative effects, and alternative="two.sided" implies selection in favor of increasingly significant effects regardless of their direction.

Beta Selection Model

When type="beta", the function can be used to fit the ‘beta selection model’ by Citkowicz and Vevea (2017). For this model, the selection function is given by \[w(p_i) = p_i^{\delta_1 - 1} \times (1 - p_i)^{\delta_2 - 1}\] where \(\delta_1 > 0\) and \(\delta_2 > 0\). The null hypothesis \(\mbox{H}_0{:}\; \delta_1 = \delta_2 = 1\) represents the case where there is no selection (at least not depending on the p-values). Note that the model allows for a non-monotonic selection function.

Half-Normal, Negative-Exponential, Logistic, and Power Selection Models

Preston et al. (2004) suggested the first three of the following selection functions:

nametypeselection function
half-normal"halfnorm"\(w(p_i) = \exp(-\delta \times p_i^2)\)
negative-exponential"negexp"\(w(p_i) = \exp(-\delta \times p_i)\)
logistic"logistic"\(w(p_i) = 2 \times \exp(-\delta \times p_i) / (1 + \exp(-\delta \times p_i))\)
power"power"\(w(p_i) = (1-p_i)^\delta\)

The power selection model is added here as it has similar properties as the models suggested by Preston et al. (2004). For all models, assume \(\delta \ge 0\), so that all functions imply a monotonically decreasing relationship between the p-value and the selection probability. For all functions, \(\mbox{H}_0{:}\; \delta = 0\) implies no selection.

Here, these functions are extended to allow for the possibility that \(w(p_i) = 1\) for p-values below a certain significance threshold denoted by \(\alpha\) (e.g., to model the case that the relative likelihood of selection is equally high for all significant studies but decreases monotonically for p-values above the significance threshold). To fit such a selection model, one should specify the \(\alpha\) value (with \(0 < \alpha < 1\)) via the steps argument. There must be at least one observed p-value below and above the chosen threshold to fit these models.

Preston et al. (2004) also suggested selection functions where the relatively likelihood of selection not only depends on the p-value, but also the precision (e.g., standard error) of the estimate (if two studies have similar p-values, it may be plausible to assume that the larger / more precise study has a higher probability of selection). These selection functions plus the corresponding power function are given by:

nametypeselection function
half-normal"halfnorm"\(w(p_i) = \exp(-\delta \times \mathrm{prec}_i \times p_i^2)\)
negative-exponential"negexp"\(w(p_i) = \exp(-\delta \times \mathrm{prec}_i \times p_i)\)
logistic"logistic"\(w(p_i) = 2 \times \exp(-\delta \times \mathrm{prec}_i \times p_i) / (1 + \exp(-\delta \times \mathrm{prec}_i \times p_i))\)
power"power"\(w(p_i) = (1-p_i)^{-\delta \times \mathrm{prec}_i}\)

where \(\mathrm{prec}_i = \sqrt{v_i}\) (i.e., the standard error of the \(i\)th study) according to Preston et al. (2004). Here, this idea is generalized to allow the user to specify the specific measure of precision to use (via the prec argument). Possible options are:

  • prec="sei" for the standard errors,

  • prec="vi" for the sampling variances,

  • prec="ninv" for the inverse of the sample sizes,

  • prec="sqrtninv" for the inverse square root of the sample sizes.

Using some function of the sample sizes as a measure of precision is only possible when information about the sample sizes is actually stored within the object passed to the selmodel function. See ‘Note’.

Note that \(\mathrm{prec}_i\) is really a measure of imprecision (with higher values corresponding to lower precision). Also, regardless of the specific measure chosen, the values are actually rescaled with \(\mathrm{prec}_i = \mathrm{prec}_i / \max(\mathrm{prec}_i)\) inside of the function, such that \(\mathrm{prec}_i = 1\) for the least precise study and \(\mathrm{prec}_i < 1\) for the remaining studies (the rescaling does not actually change the fit of the model, it only helps to improve the stability of model fitting algorithm).

One can also use the steps argument as described above in combination with these selection functions (studies with p-values below the chosen cutoff then have \(w(p_i) = 1\) regardless of their exact p-value or precision).

Negative Exponential Power Selection Model

As an extension of the half-normal and negative-exponential models, one can also choose type="negexppow" for a ‘negative exponential power selection model’. The selection function is then given by \[w(p_i) = \exp(-\delta_1 \times p_i^{1/\delta_2})\] where \(\delta_1 \ge 0\) and \(\delta_2 \ge 0\) (see Begg & Mazumdar, 1994, although here a different parameterization is used, such that increasing \(\delta_2\) leads to more severe selection).

This model affords greater flexibility in the shape of the selection function, but requires the estimation of the additional power parameter (the half-normal and negative-exponential models are therefore special cases when fixing \(\delta_2\) to 0.5 or 1, respectively). \(\mbox{H}_0{:}\; \delta_1 = 0\) again implies no selection, but so does \(\mbox{H}_0{:}\; \delta_2 = 0\).

One can again use the steps argument to specify a single significance threshold, \(\alpha\), so that \(w(p_i) = 1\) for p-values below this threshold and otherwise \(w(p_i)\) follows the selection function as given above. One can also use the prec argument to specify a measure of precision in combination with this model, which leads to the selection function \[w(p_i) = \exp(-\delta_1 \times \mathrm{prec}_i \times p_i^{1/\delta_2})\] and hence is the logical extension of the negative exponential power selection model that also incorporates some measure of precision into the selection process.

Step Function Selection Models

When type="stepfun", the function can be used to fit ‘step function models’ as described by Iyengar and Greenhouse (1988), Hedges (1992), Vevea and Hedges (1995), and Vevea and Woods (2005). For these models, one must specify one or multiple values via the steps argument, which define intervals in which the relative likelihood of selection is constant. Let \[\alpha_1 < \alpha_2 < \ldots < \alpha_c\] denote these cutpoints sorted in increasing order, with the constraint that \(\alpha_c = 1\) (if the highest value specified via steps is not 1, the function will automatically add this cutpoint), and define \(\alpha_0 = 0\). The selection function is then given by \(w(p_i) = \delta_j\) if \(\alpha_{j-1} < p_i \le \alpha_j\). To make the model identifiable, we set \(\delta_1 = 1\). The \(\delta_j\) values therefore denote the likelihood of selection in the various intervals relative to the interval for p-values between 0 and \(\alpha_1\). Hence, the null hypothesis \(\mbox{H}_0{:}\; \delta_j = 1\) for \(j = 1, \ldots, c\) implies no selection.

For example, if steps=c(.05, .10, .50, 1), then \(\delta_2\) is the likelihood of selection for p-values between .05 and .10, \(\delta_3\) is the likelihood of selection for p-values between .10 and .50, and \(\delta_4\) is the likelihood of selection for p-values between .50 and 1 relative to the likelihood of selection for p-values between 0 and .05.

There must be at least one observed p-value within each interval to fit this model. If this is not the case, an error will be issued (setting verbose=TRUE provides information about the number of p-values falling into each interval).

When applying a single cutoff in the context of a random-effects model, this model is sometimes called the ‘three-parameter selection model’ (3PSM), corresponding to the parameters \(\mu\), \(\tau^2\), and \(\delta_2\) (e.g., Carter et al., 2019; McShane et al., 2016; Pustejovsky & Rodgers, 2019). The same idea but in the context of a fixed-effects model was also described by Iyengar and Greenhouse (1988).

Note that when alternative="greater" or alternative="less" (i.e., when we assume that the relative likelihood of selection is not only related to the p-values of the studies, but also the directionality of the effects), then it would usually make sense to divide conventional levels of significance (e.g., .05) by 2 before passing these values to the steps argument. For example, if we think that studies were selected for positive effects that are significant at two-tailed \(\alpha = .05\), then we should use alternative="greater" in combination with steps=c(.025, 1).

One of the challenges when fitting this model with many cutoffs is the large number of parameters that need to be estimated (which is especially problematic when the number of studies is small). An alternative approach suggested by Vevea and Woods (2005) is to fix the \(\delta_j\) values to some a priori chosen values instead of estimating them. One can then conduct a sensitivity analysis by examining the results (e.g., the estimates of \(\mu\) and \(\tau^2\) in a random-effects model) for a variety of different sets of \(\delta_j\) values (reflecting more or less severe forms of selection). This can be done by specifying the \(\delta_j\) values via the delta argument. Table 1 in Vevea and Woods (2005) provides some illustrative examples of moderate and severe selection functions for one- and two-tailed selection. The code below creates a data frame that contains these functions.

tab <- data.frame(
   steps = c(0.005, 0.01, 0.05, 0.10, 0.25, 0.35, 0.50, 0.65, 0.75, 0.90, 0.95, 0.99, 0.995, 1),
   delta.mod.1 = c(1, 0.99, 0.95, 0.80, 0.75, 0.65, 0.60, 0.55, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50),
   delta.sev.1 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.40, 0.35, 0.30, 0.25, 0.10, 0.10, 0.10, 0.10),
   delta.mod.2 = c(1, 0.99, 0.95, 0.90, 0.80, 0.75, 0.60, 0.60, 0.75, 0.80, 0.90, 0.95, 0.99, 1.00),
   delta.sev.2 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.25, 0.25, 0.50, 0.60, 0.75, 0.90, 0.99, 1.00))

These four functions are “merely examples and should not be regarded as canonical” (Vevea & Woods, 2005).

Value

An object of class c("rma.uni","rma"). The object is a list containing the same components as a regular c("rma.uni","rma") object, but the parameter estimates are based on the selection model. Most importantly, the following elements are modified based on the selection model:

beta

estimated coefficients of the model.

se

standard errors of the coefficients.

zval

test statistics of the coefficients.

pval

corresponding p-values.

ci.lb

lower bound of the confidence intervals for the coefficients.

ci.ub

upper bound of the confidence intervals for the coefficients.

vb

variance-covariance matrix of the estimated coefficients.

tau2

estimated amount of (residual) heterogeneity. Always 0 when method="FE".

se.tau2

standard error of the estimated amount of (residual) heterogeneity.

In addition, the object contains the following additional elements:
delta

estimated selection model parameter(s).

se.delta

corresponding standard error(s).

zval.delta

corresponding test statistic(s).

pval.delta

corresponding p-value(s).

ci.lb.delta

lower bound of the confidence intervals for the parameter(s).

ci.ub.delta

upper bound of the confidence intervals for the parameter(s).

LRT

test statistic of the likelihood ratio test for the selection model parameter(s).

LRTdf

degrees of freedom for the likelihood ratio test.

LRTp

p-value for the likelihood ratio test.

LRT.tau2

test statistic of the likelihood ratio test for testing \(H_0: \tau^2 = 0\) (NA when fitting a fixed-effects model).

LRTp.tau2

p-value for the likelihood ratio test.

...

some additional elements/values.

The results of the fitted model are formatted and printed with the print.rma.uni function. The estimated selection function can be drawn with plot.rma.uni.selmodel.

Note

Model fitting is done via numerical optimization over the model parameters. By default, optim is used for the optimization. One can also chose a different optimizer via the control argument (e.g., control=list(optimizer="nlminb")). When using optim, one can set the particular method via the optmethod argument (e.g., control=list(optimizer="optim", optmethod="BFGS"), which is the default). Besides optim and nlminb, one can also choose one of the optimizers from the minqa package (i.e., uobyqa, newuoa, or bobyqa), one of the (derivative-free) algorithms from the nloptr package, the Newton-type algorithm implemented in nlm, the Hooke-Jeeves and Nelder-Mead algorithms as implemented in hjk and nmk from the dfoptim package, the quasi-Newton type optimizer ucminf from the package of the same name, or the parallelized version of the L-BFGS-B algorithm implemented in optimParallel from the package of the same name.

The optimizer name must be given as a character string (i.e., in quotes). Additional control parameters can be specified via the control argument (e.g., control=list(maxit=1000, reltol=1e-8)). For nloptr, the default is to use the BOBYQA implementation from that package with a relative convergence criterion of 1e-8 on the function value (i.e., log-likelihood), but this can be changed via the algorithm and ftop_rel arguments (e.g., control=list(optimizer="nloptr", algorithm="NLOPT_LN_SBPLX", ftol_rel=1e-6)). For optimParallel, the control argument ncpus can be used to specify the number of cores to use for the parallelization (e.g., control=list(optimizer="optimParallel", ncpus=2)). With parallel::detectCores(), one can check on the number of available cores on the local machine.

All selection models (except for type="stepfun") require repeated evaluations of an integral, which is done via adaptive quadrature as implemented in the integrate function. One can adjust the arguments of the integrate function via control element intCtrl, which is a list of named arguments (e.g., control = list(intCtrl = list(rel.tol=1e-4, subdivisions=100))).

The starting values for the fixed effects, the \(\tau^2\) value (only relevant in random/mixed-effects models), and the \(\delta\) parameter(s) are chosen automatically by the function, but one can also set the starting values manually via the control argument by specifying a vector of the appropriate length for beta.init, a single value for tau2.init, and a vector of the appropriate length for delta.init.

By default, the \(\delta\) parameter(s) are constrained to a certain range, which improves the stability of the optimization algorithm. For all models, the maximum is set to 100 and the minimum to 0 (except for type="beta", where the minimum for both parameters is 1e-05). These defaults can be changed via the control argument by specifying a vector of the appropriate length for delta.min and/or delta.max.

A difficulty with fitting the beta selection model (i.e., type="beta") is the behavior of \(w(p_i)\) when \(p_i = 0\) or \(p_i = 1\). When \(\delta_1 < 1\) or \(\delta_2 < 1\), then this leads to selection weights equal to infinity, which causes problems when computing the likelihood function. Following Citkowicz and Vevea (2017), this problem can be avoided by censoring p-values too close to 0 or 1. The specific censoring point can be set via the pval.min element of the control argument. The default for this selection model is control=list(pval.min=1e-5) (for all other selection models, the default is pval.min=0).

The variance-covariance matrix corresponding to the estimates of the fixed effects, the \(\tau^2\) value (only relevant in random/mixed-effects models), and the \(\delta\) parameter(s) is obtained by inverting the Hessian, which is numerically approximated using the hessian function. This may fail, leading to NA values for the standard errors and hence test statistics, p-values, and confidence interval bounds. One can set control argument hessianCtrl to a list of named arguments to be passed on to the method.args argument of the hessian function (the default is control=list(hessianCtrl=list(r=6))).

Information on the progress of the optimization algorithm can be obtained by setting verbose=TRUE (this won't work when using parallelization). This option is useful to determine how the model fitting is progressing. One can also set verbose to an integer (verbose=2 yields even more information and verbose=3 also show the progress visually by drawing the selection function as the optimization proceeds).

For selection functions where the prec argument is relevant, using (a function of) the sample sizes as the measure of precision (i.e., prec="ninv" or prec="sqrtninv") is only possible when information about the sample sizes is actually stored within the object passed to the selmodel function. That should automatically be the case when the observed effect sizes or outcomes were computed with the escalc function or when the observed effect sizes or outcomes were computed within the model fitting function. On the other hand, this will not be the case when rma.uni was used together with the yi and vi arguments and the yi and vi values were not computed with escalc. In that case, it is still possible to pass information about the sample sizes to the rma.uni function (e.g., use rma.uni(yi, vi, ni=ni, data=dat), where data frame dat includes a variable called ni with the sample sizes).

Finally, the automatic rescaling of the chosen precision measure can be switched off by setting scaleprec=FALSE.

Author

Wolfgang Viechtbauer wvb@metafor-project.org http://www.metafor-project.org

References

Begg, C. B., & Mazumdar, M. (1994). Operating characteristics of a rank correlation test for publication bias. Biometrics, 50(4), 1088--1101. https://doi.org/10.2307/2533446

Carter, E. C., Schönbrodt, F. D., Gervais, W. M., & Hilgard, J. (2019). Correcting for bias in psychology: A comparison of meta-analytic methods. Advances in Methods and Practices in Psychological Science, 2(2), 115--144. https://doi.org/10.1177/2515245919847196

Citkowicz, M., & Vevea, J. L. (2017). A parsimonious weight function for modeling publication bias. Psychological Methods, 22(1), 28--41. https://doi.org/10.1037/met0000119

Hedges, L. V. (1992). Modeling publication selection effects in meta-analysis. Statistical Science, 7(2), 246--255. https://doi.org/10.1214/ss/1177011364

Iyengar, S., & Greenhouse, J. B. (1988). Selection models and the file drawer problem. Statistical Science, 3(1), 109--117. https://doi.org/10.1214/ss/1177013012

McShane, B. B., Bockenholt, U., & Hansen, K. T. (2016). Adjusting for publication bias in meta-analysis: An evaluation of selection methods and some cautionary notes. Perspectives on Psychological Science, 11(5), 730--749. https://doi.org/10.1177/1745691616662243

Preston, C., Ashby, D., & Smyth, R. (2004). Adjusting for publication bias: Modelling the selection process. Journal of Evaluation in Clinical Practice, 10(2), 313--322. https://doi.org/10.1111/j.1365-2753.2003.00457.x

Pustejovsky, J. E., & Rodgers, M. A. (2019). Testing for funnel plot asymmetry of standardized mean differences. Research Synthesis Methods, 10(1), 57--71. https://doi.org/10.1002/jrsm.1332

Vevea, J. L., & Hedges, L. V. (1995). A general linear model for estimating effect size in the presence of publication bias. Psychometrika, 60(3), 419--435. https://doi.org/10.1007/BF02294384

Vevea, J. L., & Woods, C. M. (2005). Publication bias in research synthesis: Sensitivity analysis using a priori weight functions. Psychological Methods, 10(4), 428--443. https://doi.org/10.1037/1082-989X.10.4.428

See also

Examples

############################################################################ ### example from Citkowicz and Vevea (2017) for beta selection model # copy data into 'dat' and examine data dat <- dat.baskerville2012 dat
#> author year score design alloconc blind itt fumonths retention #> 1 Kottke et al. 1992 6 cct 0 1 1 19 83.0 #> 2 McBride et al. 2000 6 rct 0 0 0 18 100.0 #> 3 Stange et al. 2000 6 rct 0 0 0 24 NA #> 4 Lobo et al. 2004 6 rct 1 0 0 21 57.0 #> 5 Roetzhiem et al. 2005 6 crct 0 1 0 24 100.0 #> 6 Hogg et al. 2008 6 cct 0 0 1 6 87.0 #> 7 Aspy et al. 2008 6 cct 0 1 0 18 89.0 #> 8 Jaen et al. 2010 6 rct 0 1 0 26 86.0 #> 9 Cockburn et al. 1992 7 rct 0 0 0 3 79.0 #> 10 Modell et al. 1998 7 rct 0 0 0 12 100.0 #> 11 Engels et al. 2006 7 rct 1 0 1 12 92.0 #> 12 Aspy et al. 2008 7 rct 0 1 0 9 100.0 #> 13 Deitrich et al. 1992 18 rct 0 1 0 12 96.0 #> 14 Lobo et al. 2002 8 rct 1 0 1 21 100.0 #> 15 Bryce et al. 1995 9 rct 1 1 1 12 93.3 #> 16 Kinsinger et al. 1998 9 rct 1 1 0 18 94.0 #> 17 Solberg et al. 1998 9 rct 1 0 1 22 100.0 #> 18 Lemelin et al. 2001 9 rct 1 1 0 18 98.0 #> 19 Frijling et al. 2002 9 crct 1 1 1 21 95.0 #> 20 Frijling et al. 2003 9 crct 1 1 1 21 95.0 #> 21 Margolis et al. 2004 10 rct 1 1 1 30 100.0 #> 22 Mold et al. 2008 10 rct 1 1 1 6 100.0 #> 23 Hogg et al. 2008 12 rct 1 1 1 13 100.0 #> country outcomes duration pperf meetings hours tailor smd se #> 1 US 2 18 5.5 30.0 1.00 1 1.01 0.52 #> 2 US 4 12 NA 5.0 1.00 1 0.82 0.46 #> 3 US 35 12 20.0 4.0 1.50 1 0.59 0.23 #> 4 NL 16 21 5.0 15.0 1.00 1 0.44 0.18 #> 5 US 3 24 4.0 4.0 1.00 0 0.84 0.29 #> 6 Can 26 12 11.0 12.0 1.50 1 0.73 0.29 #> 7 US 4 18 3.0 3.0 6.00 1 1.12 0.36 #> 8 US 11 26 6.0 4.5 6.00 1 0.04 0.37 #> 9 Aus 2 2 40.0 2.0 0.25 0 0.24 0.15 #> 10 UK 1 12 13.0 3.0 1.00 0 0.32 0.40 #> 11 NL 7 5 NA 5.0 1.00 1 1.04 0.32 #> 12 US 1 9 6.0 18.0 6.00 1 1.31 0.57 #> 13 NL 16 21 20.0 15.0 1.00 1 0.59 0.29 #> 14 US 10 3 8.0 3.0 1.00 1 0.66 0.19 #> 15 US 4 12 12.0 1.0 15.00 0 0.62 0.31 #> 16 US 5 12 13.0 10.0 0.75 1 0.47 0.27 #> 17 US 10 22 11.0 4.0 3.00 1 1.08 0.32 #> 18 Can 13 18 8.0 33.0 1.75 1 0.98 0.32 #> 19 NL 7 21 20.0 15.0 1.00 0 0.26 0.18 #> 20 NL 12 21 20.0 15.0 1.00 1 0.39 0.18 #> 21 US 4 24 11.0 9.0 1.00 1 0.60 0.31 #> 22 US 6 6 8.0 18.0 4.00 1 0.94 0.53 #> 23 Can 53 12 14.0 9.0 0.75 1 0.11 0.27
# fit random-effects model res <- rma(smd, se^2, data=dat, method="ML", digits=3) res
#> #> Random-Effects Model (k = 23; tau^2 estimator: ML) #> #> tau^2 (estimated amount of total heterogeneity): 0.016 (SE = 0.024) #> tau (square root of estimated tau^2 value): 0.127 #> I^2 (total heterogeneity / total variability): 18.64% #> H^2 (total variability / sampling variability): 1.23 #> #> Test for Heterogeneity: #> Q(df = 22) = 27.552, p-val = 0.191 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> 0.555 0.063 8.775 <.001 0.431 0.679 *** #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
# funnel plot funnel(res, ylim=c(0,0.6), xlab="Standardized Mean Difference")
# fit beta selection model # \dontrun{ sel <- selmodel(res, type="beta") sel
#> #> Random-Effects Model (k = 23; tau^2 estimator: ML) #> #> tau^2 (estimated amount of total heterogeneity): 0.000 (SE = 0.001) #> tau (square root of estimated tau^2 value): 0.002 #> #> Test for Heterogeneity: #> LRT(df = 1) = 0.009, p-val = 0.925 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> 0.115 0.166 0.689 0.491 -0.212 0.441 #> #> Test for Selection Model Parameters: #> LRT(df = 2) = 7.847, p-val = 0.020 #> #> Selection Model Results: #> #> estimate se zval pval ci.lb ci.ub #> delta.1 0.473 0.235 -2.240 0.025 0.012 0.934 * #> delta.2 4.461 2.185 1.584 0.113 0.180 8.743 #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
# plot selection function plot(sel, ylim=c(0,40))# }
# fit mixed-effects meta-regression model with 'blind' dummy variable as moderator res <- rma(smd, se^2, data=dat, mods = ~ blind, method="ML", digits=3) res
#> #> Mixed-Effects Model (k = 23; tau^2 estimator: ML) #> #> tau^2 (estimated amount of residual heterogeneity): 0.016 (SE = 0.024) #> tau (square root of estimated tau^2 value): 0.128 #> I^2 (residual heterogeneity / unaccounted variability): 18.26% #> H^2 (unaccounted variability / sampling variability): 1.22 #> R^2 (amount of heterogeneity accounted for): 0.00% #> #> Test for Residual Heterogeneity: #> QE(df = 21) = 27.503, p-val = 0.155 #> #> Test of Moderators (coefficient 2): #> QM(df = 1) = 0.069, p-val = 0.792 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> intrcpt 0.573 0.093 6.192 <.001 0.392 0.754 *** #> blind -0.033 0.127 -0.263 0.792 -0.282 0.215 #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
# predicted average effect for studies that do not and that do use blinding predict(res, newmods=c(0,1))
#> #> pred se ci.lb ci.ub pi.lb pi.ub #> 1 0.573 0.093 0.392 0.754 0.264 0.882 #> 2 0.540 0.087 0.370 0.710 0.237 0.842 #>
# fit beta selection model # \dontrun{ sel <- selmodel(res, type="beta") sel
#> #> Mixed-Effects Model (k = 23; tau^2 estimator: ML) #> #> tau^2 (estimated amount of residual heterogeneity): 0.000 (SE = 0.000) #> tau (square root of estimated tau^2 value): 0.001 #> #> Test for Residual Heterogeneity: #> LRT(df = 1) = 0.000, p-val = 1.000 #> #> Test of Moderators (coefficient 2): #> QM(df = 1) = 1.201, p-val = 0.273 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> intrcpt 0.134 0.171 0.787 0.432 -0.200 0.469 #> blind -0.136 0.124 -1.096 0.273 -0.380 0.108 #> #> Test for Selection Model Parameters: #> LRT(df = 2) = 9.044, p-val = 0.011 #> #> Selection Model Results: #> #> estimate se zval pval ci.lb ci.ub #> delta.1 0.420 0.239 -2.425 0.015 0.000 0.889 * #> delta.2 5.096 2.411 1.699 0.089 0.371 9.821 . #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
predict(sel, newmods=c(0,1))# }
#> #> pred se ci.lb ci.ub pi.lb pi.ub #> 1 0.134 0.171 -0.200 0.469 -0.200 0.469 #> 2 -0.002 0.199 -0.392 0.388 -0.392 0.388 #>
############################################################################ ### example from Preston et al. (2004) # copy data into 'dat' and examine data dat <- dat.hahn2001 dat
#> study ai n1i ci n2i #> 1 Banhladesh 1995a 4 19 5 19 #> 2 Banhladesh 1996a 0 18 0 18 #> 3 CHOICE 2001 34 341 50 334 #> 4 Colombia 2000 7 71 16 69 #> 5 Egypt 1886a 6 45 5 44 #> 6 Egypt 1996b 1 94 8 96 #> 7 India 1984a 0 22 0 22 #> 8 India 2000b 11 88 12 82 #> 9 Mexico 1990a 2 82 7 84 #> 10 Panama 1982 0 33 0 30 #> 11 USA 1982 0 15 1 20 #> 12 WHO 1995 33 221 43 218
### meta-analysis of (log) odds rations using the Mantel-Haenszel method res <- rma.mh(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, digits=2, slab=study)
#> Warning: Some yi/vi values are NA.
res
#> #> Fixed-Effects Model (k = 12) #> #> I^2 (total heterogeneity / total variability): 0.00% #> H^2 (total variability / sampling variability): 0.82 #> #> Test for Heterogeneity: #> Q(df = 8) = 6.53, p-val = 0.59 #> #> Model Results (log scale): #> #> estimate se zval pval ci.lb ci.ub #> -0.49 0.14 -3.51 <.01 -0.77 -0.22 #> #> Model Results (OR scale): #> #> estimate ci.lb ci.ub #> 0.61 0.46 0.80 #> #> Cochran-Mantel-Haenszel Test: CMH = 12.00, df = 1, p-val < 0.01 #> Tarone's Test for Heterogeneity: X^2 = 7.58, df = 8, p-val = 0.48 #>
# calculate log odds ratios and corresponding sampling variances dat <- escalc(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, drop00=TRUE) dat
#> study ai n1i ci n2i yi vi #> 1 Banhladesh 1995a 4 19 5 19 -0.2921 0.5881 #> 2 Banhladesh 1996a 0 18 0 18 NA NA #> 3 CHOICE 2001 34 341 50 334 -0.4635 0.0562 #> 4 Colombia 2000 7 71 16 69 -1.0153 0.2399 #> 5 Egypt 1886a 6 45 5 44 0.1823 0.4179 #> 6 Egypt 1996b 1 94 8 96 -2.1347 1.1471 #> 7 India 1984a 0 22 0 22 NA NA #> 8 India 2000b 11 88 12 82 -0.1823 0.2015 #> 9 Mexico 1990a 2 82 7 84 -1.2910 0.6683 #> 10 Panama 1982 0 33 0 30 NA NA #> 11 USA 1982 0 15 1 20 -0.8690 2.7825 #> 12 WHO 1995 33 221 43 218 -0.3363 0.0646
# fit fixed-effects model res <- rma(yi, vi, data=dat, method="FE")
#> Warning: Studies with NAs omitted from model fitting.
# predicted odds ratio (with 95% CI) predict(res, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.63 0.48 0.83 #>
# funnel plot funnel(res, atransf=exp, at=log(c(0.01, 0.1, 1, 10, 100)), ylim=c(0,2))
# fit half-normal, negative-exponential, logistic, and power selection models sel1 <- selmodel(res, type="halfnorm", alternative="less") sel2 <- selmodel(res, type="negexp", alternative="less") sel3 <- selmodel(res, type="logistic", alternative="less") sel4 <- selmodel(res, type="power", alternative="less") # plot selection functions plot(sel1)
plot(sel2, add=TRUE, col="blue")
plot(sel3, add=TRUE, col="red")
plot(sel4, add=TRUE, col="green")
# show estimates of delta (and corresponding SEs) tab <- data.frame(delta = c(sel1$delta, sel2$delta, sel3$delta, sel4$delta), se = c(sel1$se.delta, sel2$se.delta, sel3$se.delta, sel4$se.delta)) rownames(tab) <- c("Half-normal", "Negative-exponential", "Logistic", "Power") round(tab, 2)
#> delta se #> Half-normal 3.16 2.99 #> Negative-exponential 2.66 2.35 #> Logistic 3.34 2.39 #> Power 1.46 1.39
# predicted odds ratios (with 95% CI) predict(res, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.63 0.48 0.83 #>
predict(sel1, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.71 0.49 1.04 #>
predict(sel2, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.76 0.48 1.20 #>
predict(sel3, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.74 0.49 1.13 #>
predict(sel4, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.74 0.49 1.12 #>
# fit selection models including standard error as precision measure sel1 <- selmodel(res, type="halfnorm", prec="sei", alternative="less", scaleprec=FALSE) sel2 <- selmodel(res, type="negexp", prec="sei", alternative="less", scaleprec=FALSE) sel3 <- selmodel(res, type="logistic", prec="sei", alternative="less", scaleprec=FALSE) sel4 <- selmodel(res, type="power", prec="sei", alternative="less", scaleprec=FALSE) # show estimates of delta (and corresponding SEs) tab <- data.frame(delta = c(sel1$delta, sel2$delta, sel3$delta, sel4$delta), se = c(sel1$se.delta, sel2$se.delta, sel3$se.delta, sel4$se.delta)) rownames(tab) <- c("Half-normal", "Negative-exponential", "Logistic", "Power") round(tab, 2)
#> delta se #> Half-normal 3.51 3.39 #> Negative-exponential 2.28 2.13 #> Logistic 3.02 2.32 #> Power 1.44 1.38
# predicted odds ratio (with 95% CI) predict(res, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.63 0.48 0.83 #>
predict(sel1, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.68 0.49 0.95 #>
predict(sel2, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.69 0.50 0.95 #>
predict(sel3, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.68 0.50 0.93 #>
predict(sel4, transf=exp, digits=2)
#> #> pred ci.lb ci.ub #> 0.69 0.50 0.96 #>
############################################################################ ### meta-analysis on the effect of environmental tobacco smoke on lung cancer risk # copy data into 'dat' and examine data dat <- dat.hackshaw1998 dat
#> study author year country design cases or or.lb or.ub #> 1 1 Garfinkel 1981 USA cohort 153 1.18 0.90 1.54 #> 2 2 Hirayama 1984 Japan cohort 200 1.45 1.02 2.08 #> 3 3 Butler 1988 USA cohort 8 2.02 0.48 8.56 #> 4 4 Cardenas 1997 USA cohort 150 1.20 0.80 1.60 #> 5 5 Chan 1982 Hong Kong case-control 84 0.75 0.43 1.30 #> 6 6 Correa 1983 USA case-control 22 2.07 0.81 5.25 #> 7 7 Trichopolous 1983 Greece case-control 62 2.13 1.19 3.83 #> 8 8 Buffler 1984 USA case-control 41 0.80 0.34 1.90 #> 9 9 Kabat 1984 USA case-control 24 0.79 0.25 2.45 #> 10 10 Lam 1985 Hong Kong case-control 60 2.01 1.09 3.72 #> 11 11 Garfinkel 1985 USA case-control 134 1.23 0.81 1.87 #> 12 12 Wu 1985 USA case-control 29 1.20 0.50 3.30 #> 13 13 Akiba 1986 Japan case-control 94 1.52 0.87 2.63 #> 14 14 Lee 1986 UK case-control 32 1.03 0.41 2.55 #> 15 15 Koo 1987 Hong Kong case-control 86 1.55 0.90 2.67 #> 16 16 Pershagen 1987 Sweden case-control 70 1.03 0.61 1.74 #> 17 17 Humble 1987 USA case-control 20 2.34 0.81 6.75 #> 18 18 Lam 1987 Hong Kong case-control 199 1.65 1.16 2.35 #> 19 19 Gao 1987 China case-control 246 1.19 0.82 1.73 #> 20 20 Brownson 1987 USA case-control 19 1.52 0.39 5.96 #> 21 21 Geng 1988 China case-control 54 2.16 1.08 4.29 #> 22 22 Shimizu 1988 Japan case-control 90 1.08 0.64 1.82 #> 23 23 Inoue 1988 Japan case-control 22 2.55 0.74 8.78 #> 24 24 Kalandidi 1990 Greece case-control 90 1.62 0.90 2.91 #> 25 25 Sobue 1990 Japan case-control 144 1.06 0.74 1.52 #> 26 26 Wu-Williams 1990 China case-control 417 0.79 0.62 1.02 #> 27 27 Liu 1991 China case-control 54 0.74 0.32 1.69 #> 28 28 Jockel 1991 Germany case-control 23 2.27 0.75 6.82 #> 29 29 Brownson 1992 USA case-control 431 0.97 0.78 1.21 #> 30 30 Stockwell 1992 USA case-control 210 1.60 0.80 3.00 #> 31 31 Du 1993 China case-control 75 1.19 0.66 2.13 #> 32 32 Liu 1993 China case-control 38 1.66 0.73 3.78 #> 33 33 Fontham 1994 USA case-control 651 1.26 1.04 1.54 #> 34 34 Kabat 1995 USA case-control 67 1.10 0.62 1.96 #> 35 35 Zaridze 1995 Russia case-control 162 1.66 1.12 2.45 #> 36 36 Sun 1996 China case-control 230 1.16 0.80 1.69 #> 37 37 Wang 1996 China case-control 135 1.11 0.67 1.84 #> yi vi #> 1 0.1655 0.0188 #> 2 0.3716 0.0330 #> 3 0.7031 0.5402 #> 4 0.1823 0.0313 #> 5 -0.2877 0.0797 #> 6 0.7275 0.2273 #> 7 0.7561 0.0889 #> 8 -0.2231 0.1927 #> 9 -0.2357 0.3390 #> 10 0.6981 0.0981 #> 11 0.2070 0.0456 #> 12 0.1823 0.2317 #> 13 0.4187 0.0796 #> 14 0.0296 0.2174 #> 15 0.4383 0.0770 #> 16 0.0296 0.0715 #> 17 0.8502 0.2926 #> 18 0.5008 0.0324 #> 19 0.1740 0.0363 #> 20 0.4187 0.4839 #> 21 0.7701 0.1238 #> 22 0.0770 0.0711 #> 23 0.9361 0.3982 #> 24 0.4824 0.0896 #> 25 0.0583 0.0337 #> 26 -0.2357 0.0161 #> 27 -0.3011 0.1802 #> 28 0.8198 0.3171 #> 29 -0.0305 0.0125 #> 30 0.4700 0.1137 #> 31 0.1740 0.0893 #> 32 0.5068 0.1760 #> 33 0.2311 0.0100 #> 34 0.0953 0.0862 #> 35 0.5068 0.0399 #> 36 0.1484 0.0364 #> 37 0.1044 0.0664
# fit random-effects model res <- rma(yi, vi, data=dat, method="ML") res
#> #> Random-Effects Model (k = 37; tau^2 estimator: ML) #> #> tau^2 (estimated amount of total heterogeneity): 0.0204 (SE = 0.0165) #> tau (square root of estimated tau^2 value): 0.1427 #> I^2 (total heterogeneity / total variability): 27.62% #> H^2 (total variability / sampling variability): 1.38 #> #> Test for Heterogeneity: #> Q(df = 36) = 47.4979, p-val = 0.0952 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> 0.2171 0.0486 4.4712 <.0001 0.1219 0.3123 *** #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
# funnel plot funnel(res, atransf=exp, at=log(c(0.25, 0.5, 1, 2, 4, 8)), ylim=c(0,0.8))
# step function selection model # \dontrun{ sel <- selmodel(res, type="stepfun", steps=c(0.05, 0.10, 0.50, 1.00)) sel
#> #> Random-Effects Model (k = 37; tau^2 estimator: ML) #> #> tau^2 (estimated amount of total heterogeneity): 0.0307 (SE = 0.0271) #> tau (square root of estimated tau^2 value): 0.1753 #> #> Test for Heterogeneity: #> LRT(df = 1) = 3.9537, p-val = 0.0468 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> 0.1217 0.1299 0.9373 0.3486 -0.1328 0.3763 #> #> Test for Selection Model Parameters: #> LRT(df = 3) = 7.0661, p-val = 0.0698 #> #> Selection Model Results: #> #> k estimate se zval pval ci.lb ci.ub #> 0 < p <= 0.05 7 1.0000 --- --- --- --- --- #> 0.05 < p <= 0.1 8 2.4221 1.6609 0.8562 0.3919 0.0000 5.6773 #> 0.1 < p <= 0.5 16 0.9775 0.8204 -0.0274 0.9782 0.0000 2.5855 #> 0.5 < p <= 1 6 0.3967 0.4692 -1.2857 0.1986 0.0000 1.3164 #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
# plot selection function plot(sel)# }
############################################################################ ### validity of student ratings example from Vevea & Woods (2005) # copy data into 'dat' and examine data dat <- dat.cohen1981 dat
#> study sample #> 1 Bolton et al. 1979 General psychology #> 2 Bryson 1974 College algebra #> 3 Centra 1977 General biology #> 4 Centra 1977 General psychology #> 5 Crooks & Smock 1974 General physics #> 6 Doyle & Crichton 1978 Introductory communications #> 7 Doyle & Whitely 1974 Beginning French #> 8 Elliott 1950 General chemistry #> 9 Ellis & Rickard 1977 General psychology #> 10 Frey et al. 1975 Introductory calculus #> 11 Greenwood et al. 1976 Analytic geometry and calculus #> 12 Hoffman 1978 Introductory math #> 13 McKeachie et al. 1971 General psychology #> 14 Morsh et al. 1956 Aircraft mechanics #> 15 Remmers et al. 1949 General chemistry #> 16 Sullivan & Skanes 1974 First-year science #> 17 Sullivan & Skanes 1974 Introductory psychology #> 18 Sullivan & Skanes 1974 First-year math #> 19 Sullivan & Skanes 1974 First-year biology #> 20 Wherry 1952 Introductory psychology #> control ni ri #> 1 Random assignment 10 0.68 #> 2 CIAT pretest 20 0.56 #> 3 Final exam for prior semester of same course 13 0.23 #> 4 Final exam for prior semester of same course 22 0.64 #> 5 Math pretest 28 0.49 #> 6 PSAT verbal score 12 -0.04 #> 7 Minnesota Scholastic Aptitude Test 12 0.49 #> 8 College entrance exam scores 36 0.33 #> 9 Random assignment 19 0.58 #> 10 SAT math scores 12 0.18 #> 11 Achievement pretest 36 -0.11 #> 12 SAT math scores 75 0.27 #> 13 Intelligence test scores 33 0.26 #> 14 Hydraulics pretest 121 0.40 #> 15 Orientation test scores 37 0.49 #> 16 Random assignment 14 0.51 #> 17 Random assignment 40 0.40 #> 18 Random assignment 16 0.34 #> 19 Random assignment 14 0.42 #> 20 Ability scores 20 0.16
# calculate r-to-z transformed correlations and corresponding sampling variances dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat[c(1,4,5)]) dat
#> study ni ri yi vi #> 1 Bolton et al. 1979 10 0.68 0.8291 0.1429 #> 2 Bryson 1974 20 0.56 0.6328 0.0588 #> 3 Centra 1977 13 0.23 0.2342 0.1000 #> 4 Centra 1977 22 0.64 0.7582 0.0526 #> 5 Crooks & Smock 1974 28 0.49 0.5361 0.0400 #> 6 Doyle & Crichton 1978 12 -0.04 -0.0400 0.1111 #> 7 Doyle & Whitely 1974 12 0.49 0.5361 0.1111 #> 8 Elliott 1950 36 0.33 0.3428 0.0303 #> 9 Ellis & Rickard 1977 19 0.58 0.6625 0.0625 #> 10 Frey et al. 1975 12 0.18 0.1820 0.1111 #> 11 Greenwood et al. 1976 36 -0.11 -0.1104 0.0303 #> 12 Hoffman 1978 75 0.27 0.2769 0.0139 #> 13 McKeachie et al. 1971 33 0.26 0.2661 0.0333 #> 14 Morsh et al. 1956 121 0.40 0.4236 0.0085 #> 15 Remmers et al. 1949 37 0.49 0.5361 0.0294 #> 16 Sullivan & Skanes 1974 14 0.51 0.5627 0.0909 #> 17 Sullivan & Skanes 1974 40 0.40 0.4236 0.0270 #> 18 Sullivan & Skanes 1974 16 0.34 0.3541 0.0769 #> 19 Sullivan & Skanes 1974 14 0.42 0.4477 0.0909 #> 20 Wherry 1952 20 0.16 0.1614 0.0588
# fit random-effects model res <- rma(yi, vi, data=dat, method="ML", digits=3) res
#> #> Random-Effects Model (k = 20; tau^2 estimator: ML) #> #> tau^2 (estimated amount of total heterogeneity): 0.001 (SE = 0.009) #> tau (square root of estimated tau^2 value): 0.034 #> I^2 (total heterogeneity / total variability): 2.86% #> H^2 (total variability / sampling variability): 1.03 #> #> Test for Heterogeneity: #> Q(df = 19) = 20.974, p-val = 0.338 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> 0.380 0.045 8.505 <.001 0.292 0.468 *** #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
# predicted average correlation (with 95% CI) predict(res, transf=transf.ztor)
#> #> pred ci.lb ci.ub pi.lb pi.ub #> 0.363 0.284 0.436 0.263 0.455 #>
# funnel plot funnel(res, ylim=c(0,0.4))
# selection functions from Vevea & Woods (2005) tab <- data.frame( steps = c(0.005, 0.01, 0.05, 0.10, 0.25, 0.35, 0.50, 0.65, 0.75, 0.90, 0.95, 0.99, 0.995, 1), delta.mod.1 = c(1, 0.99, 0.95, 0.80, 0.75, 0.65, 0.60, 0.55, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50), delta.sev.1 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.40, 0.35, 0.30, 0.25, 0.10, 0.10, 0.10, 0.10), delta.mod.2 = c(1, 0.99, 0.95, 0.90, 0.80, 0.75, 0.60, 0.60, 0.75, 0.80, 0.90, 0.95, 0.99, 1.00), delta.sev.2 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.25, 0.25, 0.50, 0.60, 0.75, 0.90, 0.99, 1.00)) # apply step function model with a priori chosen selection weights # \dontrun{ sel <- lapply(tab[-1], function(delta) selmodel(res, type="stepfun", steps=tab$steps, delta=delta)) # estimates (transformed correlation) and tau^2 values sav <- data.frame(estimate = round(c(res$beta, sapply(sel, function(x) x$beta)), 2), varcomp = round(c(res$tau2, sapply(sel, function(x) x$tau2)), 3)) sav# }
#> estimate varcomp #> 0.38 0.001 #> delta.mod.1 0.35 0.005 #> delta.sev.1 0.32 0.010 #> delta.mod.2 0.36 0.003 #> delta.sev.2 0.33 0.006
############################################################################