Function to fit regression models based on correlation and covariance matrices.

matreg(y, x, R, n, V, cov=FALSE, means, ztor=FALSE,
       nearpd=FALSE, level=95, digits, ...)

Arguments

y

index (or name given as a character string) of the outcome variable.

x

indices (or names given as a character vector) of the predictor variables.

R

correlation or covariance matrix (or only the lower triangular part including the diagonal).

n

sample size based on which the elements in the correlation/covariance matrix were computed.

V

variance-covariance matrix of the lower triangular elements of the correlation/covariance matrix. Either V or n should be specified, not both. See ‘Details’.

cov

logical to specify whether R is a covariance matrix (the default is FALSE).

means

optional vector to specify the means of the variables (only relevant when cov=TRUE).

ztor

logical to specify whether R is a matrix of r-to-z transformed correlations and hence should be back-transformed to raw correlations (the default is FALSE). See ‘Details’.

nearpd

logical to specify whether the nearPD function from the Matrix package should be used when the \(R_{x,x}\) matrix cannot be inverted. See ‘Note’.

level

numeric value between 0 and 100 to specify the confidence interval level (the default is 95; see here for details).

digits

optional integer to specify the number of decimal places to which the printed results should be rounded.

...

other arguments.

Details

Let \(R\) be a \(p \times p\) correlation or covariance matrix. Let \(y\) denote the row/column of the outcome variable and \(x\) the row(s)/column(s) of the predictor variable(s) in this matrix. Let \(R_{x,x}\) and \(R_{x,y}\) denote the corresponding submatrices of \(R\). Then \[b = R_{x,x}^{-1} R_{x,y}\] yields the standardized or raw regression coefficients (depending on whether \(R\) is a correlation or covariance matrix, respectively) when regressing the outcome variable on the predictor variable(s).

The \(R\) matrix may be computed based on a single sample of \(n\) subjects. In this case, one should specify the sample size via argument n. The variance-covariance matrix of the standardized regression coefficients is then given by \(\mbox{Var}[b] = \mbox{MSE} \times R_{x,x}^{-1}\), where \(\mbox{MSE} = (1 - b'R_{x,y}) / (n - m)\) and \(m\) denotes the number of predictor variables. The standard errors are then given by the square root of the diagonal elements of \(\mbox{Var}[b]\). Test statistics (in this case, t-statistics) and the corresponding p-values can then be computed as in a regular regression analysis. When \(R\) is a covariance matrix, one should set cov=TRUE and specify the means of the \(p\) variables via argument means to obtain raw regression coefficients including the intercept and corresponding standard errors.

Alternatively, \(R\) may be the result of a meta-analysis of correlation coefficients. In this case, the elements in \(R\) are pooled correlation coefficients and the variance-covariance matrix of these pooled coefficients should be specified via argument V. The order of elements in V should correspond to the order of elements in the lower triangular part of \(R\) column-wise. For example, if \(R\) is a \(4 \times 4\) matrix of the form: \[\begin{bmatrix} 1 & & & \\ r_{21} & 1 & & \\ r_{31} & r_{32} & 1 & \\ r_{41} & r_{42} & r_{43} & 1 \end{bmatrix}\] then the elements are \(r_{21}\), \(r_{31}\), \(r_{41}\), \(r_{32}\), \(r_{42}\), and \(r_{43}\) and hence V should be a \(6 \times 6\) variance-covariance matrix of these elements in this order. The variance-covariance matrix of the standardized regression coefficients (i.e., \(\mbox{Var}[b]\)) is then computed as a function of V as described in Becker (1992) using the multivariate delta method. The standard errors are then again given by the square root of the diagonal elements of \(\mbox{Var}[b]\). Test statistics (in this case, z-statistics) and the corresponding p-values can then be computed in the usual manner.

In case \(R\) is the result of a meta-analysis of Fisher r-to-z transformed correlation coefficients (and hence V is then the corresponding variance-covariance matrix of these pooled transformed coefficients), one should set argument ztor=TRUE, so that the appropriate back-transformation is then applied to R (and V) within the function.

Finally, \(R\) may be a covariance matrix based on a meta-analysis (e.g., the estimated variance-covariance matrix of the random effects in a multivariate model). In this case, one should set cov=TRUE and V should again be the variance-covariance matrix of the elements in \(R\), but now including the diagonal. Hence, if \(R\) is a \(4 \times 4\) matrix of the form: \[\begin{bmatrix} \tau_1^2 & & & \\ \tau_{21} & \tau_2^2 & & \\ \tau_{31} & \tau_{32} & \tau_3^2 & \\ \tau_{41} & \tau_{42} & \tau_{43} & \tau_4^2 \end{bmatrix}\] then the elements are \(\tau^2_1\), \(\tau_{21}\), \(\tau_{31}\), \(\tau_{41}\), \(\tau^2_2\), \(\tau_{32}\), \(\tau_{42}\), \(\tau^2_3\), \(\tau_{43}\), and \(\tau^2_4\), and hence V should be a \(10 \times 10\) variance-covariance matrix of these elements in this order. Argument means can then again be used to specify the means of the variables.

Value

An object of class "matreg". The object is a list containing the following components:

tab

a data frame with the estimated model coefficients, standard errors, test statistics, degrees of freedom (only for t-tests), p-values, and lower/upper confidence interval bounds.

vb

the variance-covariance matrix of the estimated model coefficients.

...

some additional elements/values.

The results are formatted and printed with the print function. Extractor functions include coef and vcov.

Note

Only the lower triangular part of R (and V if it is specified) is used in the computations.

If \(R_{x,x}\) is not invertible, an error will be issued. In this case, one can set argument nearpd=TRUE, in which case the nearPD function from the Matrix package will be used to find the nearest positive semi-definite matrix, which should be invertible. The results should be treated with caution when this is done.

When \(R\) is a covariance matrix with V and means specified, the means are treated as known constants when estimating the standard error of the intercept.

References

Becker, B. J. (1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341–362. https://doi.org/10.3102/10769986017004341

Becker, B. J. (1995). Corrections to "Using results from replicated studies to estimate linear models". Journal of Educational and Behavioral Statistics, 20(1), 100–102. https://doi.org/10.3102/10769986020001100

Becker, B. J., & Aloe, A. (2019). Model-based meta-analysis and related approaches. In H. Cooper, L. V. Hedges, & J. C. Valentine (Eds.), The handbook of research synthesis and meta-analysis (3rd ed., pp. 339–363). New York: Russell Sage Foundation.

See also

rma.mv for a function to meta-analyze multiple correlation coefficients that can be used to construct an \(R\) matrix.

rcalc for a function to construct the variance-covariance matrix of dependent correlation coefficients.

Examples

############################################################################

### first an example unrelated to meta-analysis, simply demonstrating that
### one can obtain the same results from lm() and matreg()

### fit a regression model with lm() to the 'mtcars' dataset
res <- lm(mpg ~ hp + wt + am, data=mtcars)
summary(res)
#> 
#> Call:
#> lm(formula = mpg ~ hp + wt + am, data = mtcars)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.4221 -1.7924 -0.3788  1.2249  5.5317 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 34.002875   2.642659  12.867 2.82e-13 ***
#> hp          -0.037479   0.009605  -3.902 0.000546 ***
#> wt          -2.878575   0.904971  -3.181 0.003574 ** 
#> am           2.083710   1.376420   1.514 0.141268    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 2.538 on 28 degrees of freedom
#> Multiple R-squared:  0.8399,	Adjusted R-squared:  0.8227 
#> F-statistic: 48.96 on 3 and 28 DF,  p-value: 2.908e-11
#> 

### covariance matrix of the dataset
S <- cov(mtcars)

### fit the same regression model using matreg()
res <- matreg(y="mpg", x=c("hp","wt","am"), R=S, cov=TRUE,
              means=colMeans(mtcars), n=nrow(mtcars))
summary(res)
#> 
#>          estimate      se     tval  df    pval    ci.lb    ci.ub      
#> intrcpt   34.0029  2.6427  12.8669  28  <.0001  28.5896  39.4161  *** 
#> hp        -0.0375  0.0096  -3.9018  28  0.0005  -0.0572  -0.0178  *** 
#> wt        -2.8786  0.9050  -3.1808  28  0.0036  -4.7323  -1.0248   ** 
#> am         2.0837  1.3764   1.5139  28  0.1413  -0.7358   4.9032      
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 2.5375 on 28 degrees of freedom
#> Multiple R-squared: 0.8399,  Adjusted R-squared: 0.8227
#> F-statistic: 48.9600 on 3 and 28 DF,  p-value: <.0001
#> 

### copy the 'mtcars' dataset to 'dat' and standardize all variables
dat <- mtcars
dat[] <- scale(dat)

### fit a regression model with lm() to obtain standardized regression coefficients ('betas')
res <- lm(mpg ~ 0 + hp + wt + am, data=dat)
summary(res)
#> 
#> Call:
#> lm(formula = mpg ~ 0 + hp + wt + am, data = dat)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.56779 -0.29739 -0.06285  0.20324  0.91783 
#> 
#> Coefficients:
#>    Estimate Std. Error t value Pr(>|t|)    
#> hp  -0.4264     0.1074  -3.971 0.000433 ***
#> wt  -0.4673     0.1444  -3.237 0.003017 ** 
#> am   0.1725     0.1120   1.541 0.134242    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 0.4137 on 29 degrees of freedom
#> Multiple R-squared:  0.8399,	Adjusted R-squared:  0.8233 
#> F-statistic: 50.71 on 3 and 29 DF,  p-value: 1.183e-11
#> 

### correlation matrix of the dataset
R <- cor(mtcars)

### fit the same regression model using matreg()
res <- matreg(y="mpg", x=c("hp","wt","am"), R=R, n=nrow(mtcars))
summary(res)
#> 
#>     estimate      se     tval  df    pval    ci.lb    ci.ub      
#> hp   -0.4264  0.1074  -3.9709  29  0.0004  -0.6460  -0.2068  *** 
#> wt   -0.4673  0.1444  -3.2372  29  0.0030  -0.7626  -0.1721   ** 
#> am    0.1725  0.1120   1.5407  29  0.1342  -0.0565   0.4015      
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 0.4137 on 29 degrees of freedom
#> Multiple R-squared: 0.8399,  Adjusted R-squared: 0.8233
#> F-statistic: 50.7086 on 3 and 29 DF,  p-value: <.0001
#> 

### note: the standard errors of the betas should not be used to construct CIs
### as they assume that the null hypothesis (H0: beta_j = 0) is true

### construct the var-cov matrix of correlations in R
V <- rcalc(R, ni=nrow(mtcars))$V

### fit the same regression model using matreg() but now supply V
res <- matreg(y="mpg", x=c("hp","wt","am"), R=R, V=V)
summary(res)
#> 
#>     estimate      se     zval    pval    ci.lb    ci.ub      
#> hp   -0.4264  0.1060  -4.0218  <.0001  -0.6341  -0.2186  *** 
#> wt   -0.4673  0.1378  -3.3920  0.0007  -0.7374  -0.1973  *** 
#> am    0.1725  0.1097   1.5724  0.1158  -0.0425   0.3876      
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> R^2: 0.8399, QM(df = 3) = 758.2210, p-val < .0001
#> 

### the standard errors computed in this way can now be used to construct
### CIs for the betas (here, the difference is relatively small)

############################################################################

### copy data into 'dat'
dat <- dat.craft2003

### construct dataset and var-cov matrix of the correlations
tmp <- rcalc(ri ~ var1 + var2 | study, ni=ni, data=dat)
V <- tmp$V
dat <- tmp$dat

### turn var1.var2 into a factor with the desired order of levels
dat$var1.var2 <- factor(dat$var1.var2,
   levels=c("acog.perf", "asom.perf", "conf.perf", "acog.asom", "acog.conf", "asom.conf"))

### multivariate random-effects model
res <- rma.mv(yi, V, mods = ~ var1.var2 - 1, random = ~ var1.var2 | study, struct="UN", data=dat)
#> Warning: 9 rows with NAs omitted from model fitting.
res
#> 
#> Multivariate Meta-Analysis Model (k = 51; method: REML)
#> 
#> Variance Components:
#> 
#> outer factor: study     (nlvls = 9)
#> inner factor: var1.var2 (nlvls = 6)
#> 
#>             estim    sqrt  k.lvl  fixed      level 
#> tau^2.1    0.1611  0.4014      9     no  acog.perf 
#> tau^2.2    0.0604  0.2459      9     no  asom.perf 
#> tau^2.3    0.0468  0.2163      8     no  conf.perf 
#> tau^2.4    0.0047  0.0683      9     no  acog.asom 
#> tau^2.5    0.0125  0.1119      8     no  acog.conf 
#> tau^2.6    0.0111  0.1052      8     no  asom.conf 
#> 
#>            rho.acg.p  rho.asm.p  rho.cnf.  rho.acg.s  rho.acg.c  rho.asm.c    acg.p  asm.p  cnf. 
#> acog.perf          1                                                              -      9     8 
#> asom.perf     0.9497          1                                                  no      -     8 
#> conf.perf    -0.6178    -0.5969         1                                        no     no     - 
#> acog.asom     0.5491     0.4604   -0.9345          1                             no     no    no 
#> acog.conf     0.0432    -0.0495    0.7023    -0.6961          1                  no     no    no 
#> asom.conf     0.3532     0.2688   -0.1311    -0.0891     0.4193          1       no     no    no 
#>            acg.s  acg.c  asm.c 
#> acog.perf      9      8      8 
#> asom.perf      9      8      8 
#> conf.perf      8      8      8 
#> acog.asom      -      8      8 
#> acog.conf     no      -      8 
#> asom.conf     no     no      - 
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 45) = 334.8358, p-val < .0001
#> 
#> Test of Moderators (coefficients 1:6):
#> QM(df = 6) = 596.7711, p-val < .0001
#> 
#> Model Results:
#> 
#>                     estimate      se     zval    pval    ci.lb    ci.ub      
#> var1.var2acog.perf   -0.0600  0.1408  -0.4264  0.6698  -0.3359   0.2159      
#> var1.var2asom.perf   -0.1423  0.0917  -1.5527  0.1205  -0.3220   0.0373      
#> var1.var2conf.perf    0.3167  0.0847   3.7393  0.0002   0.1507   0.4827  *** 
#> var1.var2acog.asom    0.5671  0.0367  15.4640  <.0001   0.4953   0.6390  *** 
#> var1.var2acog.conf   -0.4888  0.0509  -9.6048  <.0001  -0.5886  -0.3891  *** 
#> var1.var2asom.conf   -0.4750  0.0506  -9.3901  <.0001  -0.5741  -0.3758  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### restructure estimated mean correlations into a 4x4 matrix
R <- vec2mat(coef(res))
rownames(R) <- colnames(R) <- c("perf", "acog", "asom", "conf")
round(R, digits=3)
#>        perf   acog   asom   conf
#> perf  1.000 -0.060 -0.142  0.317
#> acog -0.060  1.000  0.567 -0.489
#> asom -0.142  0.567  1.000 -0.475
#> conf  0.317 -0.489 -0.475  1.000

### check that order in vcov(res) corresponds to order in R
round(vcov(res), digits=4)
#>                    var1.var2acog.perf var1.var2asom.perf var1.var2conf.perf var1.var2acog.asom
#> var1.var2acog.perf             0.0198             0.0115            -0.0069             0.0017
#> var1.var2asom.perf             0.0115             0.0084            -0.0043             0.0009
#> var1.var2conf.perf            -0.0069            -0.0043             0.0072            -0.0017
#> var1.var2acog.asom             0.0017             0.0009            -0.0017             0.0013
#> var1.var2acog.conf             0.0004            -0.0002             0.0023            -0.0009
#> var1.var2asom.conf             0.0018             0.0010            -0.0004            -0.0004
#>                    var1.var2acog.conf var1.var2asom.conf
#> var1.var2acog.perf             0.0004             0.0018
#> var1.var2asom.perf            -0.0002             0.0010
#> var1.var2conf.perf             0.0023            -0.0004
#> var1.var2acog.asom            -0.0009            -0.0004
#> var1.var2acog.conf             0.0026             0.0011
#> var1.var2asom.conf             0.0011             0.0026

### fit regression model with 'perf' as outcome and 'acog', 'asom', and 'conf' as predictors
matreg(1, 2:4, R=R, V=vcov(res))
#> 
#>       estimate      se     zval    pval    ci.lb   ci.ub      
#> acog    0.1482  0.1566   0.9465  0.3439  -0.1587  0.4550      
#> asom   -0.0536  0.0768  -0.6979  0.4852  -0.2043  0.0970      
#> conf    0.3637  0.0910   3.9985  <.0001   0.1854  0.5419  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### can also specify variable names
matreg("perf", c("acog","asom","conf"), R=R, V=vcov(res))
#> 
#>       estimate      se     zval    pval    ci.lb   ci.ub      
#> acog    0.1482  0.1566   0.9465  0.3439  -0.1587  0.4550      
#> asom   -0.0536  0.0768  -0.6979  0.4852  -0.2043  0.0970      
#> conf    0.3637  0.0910   3.9985  <.0001   0.1854  0.5419  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### repeat the above but with r-to-z transformed correlations
dat <- dat.craft2003
tmp <- rcalc(ri ~ var1 + var2 | study, ni=ni, data=dat, rtoz=TRUE)
V <- tmp$V
dat <- tmp$dat
dat$var1.var2 <- factor(dat$var1.var2,
   levels=c("acog.perf", "asom.perf", "conf.perf", "acog.asom", "acog.conf", "asom.conf"))
res <- rma.mv(yi, V, mods = ~ var1.var2 - 1, random = ~ var1.var2 | study, struct="UN", data=dat)
#> Warning: 9 rows with NAs omitted from model fitting.
R <- vec2mat(coef(res))
rownames(R) <- colnames(R) <- c("perf", "acog", "asom", "conf")
matreg(1, 2:4, R=R, V=vcov(res), ztor=TRUE)
#> 
#>       estimate      se     zval    pval    ci.lb   ci.ub      
#> acog    0.1362  0.1697   0.8023  0.4224  -0.1965  0.4688      
#> asom   -0.0678  0.0761  -0.8900  0.3735  -0.2170  0.0814      
#> conf    0.3665  0.0934   3.9248  <.0001   0.1835  0.5496  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

############################################################################

### a different example based on van Houwelingen et al. (2002)

### create dataset in long format
dat.long <- to.long(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg,
                    data=dat.colditz1994, append=FALSE)
dat.long <- escalc(measure="PLO", xi=out1, mi=out2, data=dat.long)
levels(dat.long$group) <- c("CON", "EXP")
dat.long
#> 
#>    study group out1  out2      yi     vi 
#> 1      1   EXP    4   119 -3.3928 0.2584 
#> 2      1   CON   11   128 -2.4541 0.0987 
#> 3      2   EXP    6   300 -3.9120 0.1700 
#> 4      2   CON   29   274 -2.2458 0.0381 
#> 5      3   EXP    3   228 -4.3307 0.3377 
#> 6      3   CON   11   209 -2.9444 0.0957 
#> 7      4   EXP   62 13536 -5.3860 0.0162 
#> 8      4   CON  248 12619 -3.9295 0.0041 
#> 9      5   EXP   33  5036 -5.0279 0.0305 
#> 10     5   CON   47  5761 -4.8087 0.0215 
#> 11     6   EXP  180  1361 -2.0230 0.0063 
#> 12     6   CON  372  1079 -1.0649 0.0036 
#> 13     7   EXP    8  2537 -5.7593 0.1254 
#> 14     7   CON   10   619 -4.1255 0.1016 
#> 15     8   EXP  505 87886 -5.1592 0.0020 
#> 16     8   CON  499 87892 -5.1713 0.0020 
#> 17     9   EXP   29  7470 -5.5514 0.0346 
#> 18     9   CON   45  7232 -5.0796 0.0224 
#> 19    10   EXP   17  1699 -4.6046 0.0594 
#> 20    10   CON   65  1600 -3.2034 0.0160 
#> 21    11   EXP  186 50448 -5.6030 0.0054 
#> 22    11   CON  141 27197 -5.2621 0.0071 
#> 23    12   EXP    5  2493 -6.2118 0.2004 
#> 24    12   CON    3  2338 -6.6584 0.3338 
#> 25    13   EXP   27 16886 -6.4384 0.0371 
#> 26    13   CON   29 17825 -6.4211 0.0345 
#> 

### fit bivariate model
res <- rma.mv(yi, vi, mods = ~ group - 1, random = ~ group | study, struct="UN",
              data=dat.long, method="ML")
res
#> 
#> Multivariate Meta-Analysis Model (k = 26; method: ML)
#> 
#> Variance Components:
#> 
#> outer factor: study (nlvls = 13)
#> inner factor: group (nlvls = 2)
#> 
#>             estim    sqrt  k.lvl  fixed  level 
#> tau^2.1    2.4073  1.5516     13     no    CON 
#> tau^2.2    1.4314  1.1964     13     no    EXP 
#> 
#>      rho.CON  rho.EXP    CON  EXP 
#> CON        1               -   13 
#> EXP   0.9467        1     no    - 
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 24) = 5270.3863, p-val < .0001
#> 
#> Test of Moderators (coefficients 1:2):
#> QM(df = 2) = 292.4633, p-val < .0001
#> 
#> Model Results:
#> 
#>           estimate      se      zval    pval    ci.lb    ci.ub      
#> groupCON   -4.0960  0.4347   -9.4226  <.0001  -4.9480  -3.2440  *** 
#> groupEXP   -4.8337  0.3396  -14.2329  <.0001  -5.4994  -4.1681  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### regression of log(odds)_EXP on log(odds)_CON
matreg(y=2, x=1, R=res$G, cov=TRUE, means=coef(res), n=res$g.levels.comb.k)
#> 
#>          estimate      se     tval  df    pval    ci.lb    ci.ub      
#> intrcpt   -1.8437  0.3265  -5.6477  11  0.0001  -2.5623  -1.1252  *** 
#> CON        0.7300  0.0749   9.7467  11  <.0001   0.5651   0.8948  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### but the SE of the CON coefficient is not computed correctly, since above we treat res$G as if
### it was a var-cov matrix computed from raw data based on res$g.levels.comb.k (= 13) data points

### fit bivariate model and get the var-cov matrix of the estimates in res$G
res <- rma.mv(yi, vi, mods = ~ group - 1, random = ~ group | study, struct="UN",
              data=dat.long, method="ML", cvvc="varcov", control=list(nearpd=TRUE))

### now use res$vvc as the var-cov matrix of the estimates in res$G
matreg(y=2, x=1, R=res$G, cov=TRUE, means=coef(res), V=res$vvc)
#> 
#>          estimate      se     zval    pval    ci.lb    ci.ub      
#> intrcpt   -1.8437  0.3548  -5.1967  <.0001  -2.5391  -1.1484  *** 
#> CON        0.7300  0.0866   8.4276  <.0001   0.5602   0.8998  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>