Cover of the book

Chapter 3: Computing Effect Sizes for Meta-Analysis

Repeated measures

###################################################
### code chunk number 2: RM1
###################################################
## Library for the SEM
library("lavaan")

## Sample covariance matrix on pre- and post-test scores
lower <- '10
          8 12'
( Cov <- getCov(lower, diag=TRUE, names=c("x_pre","x_post")) )
##        x_pre x_post
## x_pre     10      8
## x_post     8     12
## Sample means for the pre- and post-test scores
Mean <- c(10, 13)

## Sample size
N <- 50

model1 <- '# Label the sds with sd_pre and sd_post
           eta_pre =~ sd_pre*x_pre
           eta_post =~ sd_post*x_post
           # r: correlation betwen pre- and post-test
           eta_pre ~~ r*eta_post
           # Fix the error variances at 0
           x_pre ~~ 0*x_pre
           x_post ~~ 0*x_post
           # Label the means with m_pre and m_post 
           x_pre ~ m_pre*1
           x_post ~ m_post*1
           # Calculate the effect sizes
           # Diff: change score
           Diff := m_post - m_pre
           # SMD.cs sd on change score as the standardizer
           SMD.cs := (m_post - m_pre)/sqrt(sd_pre^2+sd_post^2
                      -2*sd_pre*sd_post*r)
           # SMD.pre: sd_pre as the standardizer
           SMD.pre := (m_post - m_pre)/sd_pre'

## Fit the model
fit1 <- cfa(model1, sample.cov=Cov, sample.mean=Mean, 
            sample.nobs=N, std.lv=TRUE, 
            sample.cov.rescale=FALSE)

## Display the summary
## summary(fit1)

## Display the selected output
parameterEstimates(fit1)[c(12,13,14), -c(1,2,3)]
##      label   est    se     z pvalue ci.lower ci.upper
## 12    Diff 3.000 0.346 8.660      0    2.321    3.679
## 13  SMD.cs 1.225 0.187 6.547      0    0.858    1.591
## 14 SMD.pre 0.949 0.145 6.547      0    0.665    1.233
###################################################
### code chunk number 4: RM3
###################################################
model2 <- '# Label the common sd
           eta_pre =~ sd*x_pre
           eta_post =~ sd*x_post
           # r: correlation betwen pre- and post-test
           eta_pre ~~ r*eta_post
           # Fix the error variances at 0
           x_pre ~~ 0*x_pre
           x_post ~~ 0*x_post
           # Label the means with m_pre and m_post 
           x_pre ~ m_pre*1
           x_post ~ m_post*1
           # Calculate the effect sizes
           # Common sd
           SMD.common := (m_post-m_pre)/sd'

## Fit the model
fit2 <- cfa(model2, sample.cov=Cov, sample.mean=Mean, 
            sample.nobs=N, std.lv=TRUE, 
            sample.cov.rescale=FALSE)

## Display the selected output
parameterEstimates(fit2)[12, -c(1,2,3)]
##         label   est    se     z pvalue ci.lower ci.upper
## 12 SMD.common 0.905 0.131 6.904      0    0.648    1.161

Multiple treatment studies

###################################################
### code chunk number 6: MT1
###################################################
## Group 1 (control group): variance 
var1 <- matrix(10, dimnames=list("x","x"))
## Group 2 (treatment 1): variance 
var2 <- matrix(11, dimnames=list("x","x"))
## Group 3 (treatment 2): variance 
var3 <- matrix(12, dimnames=list("x","x"))

## Convert variances into a list
Var <- list(var1, var2, var3)

## Means for the groups
Mean <- list(5, 7, 9)

## Sample sizes for the groups
N <- c(50, 52, 53)

## Assuming homogeneity of variances by using the same label "s2"
model3 <- 'x ~~ c("s2", "s2", "s2")*x
           x ~ c("m1", "m2", "m3")*1
           # SMD for treatment 1
           MT1 := (m2-m1)/sqrt(s2)
           # SMD for treatment 2
           MT2 := (m3-m1)/sqrt(s2)'

fit3 <- sem(model3, sample.cov=Var, sample.mean=Mean, 
            sample.nobs=N, sample.cov.rescale=FALSE)

## Obtain the free parameters in the model
( x <- fit3@Fit@x )
## [1] 11.01935  5.00000 11.01935  7.00000 11.01935  9.00000
## Obtain the sampling covariance matrix of the parameter estimates
( VCOV <- vcov(fit3) )
##    s2    m1    s2    m2    s2    m3   
## s2 1.567                              
## m1 0.000 0.220                        
## s2 1.567 0.000 1.567                  
## m2 0.000 0.000 0.000 0.212            
## s2 1.567 0.000 1.567 0.000 1.567      
## m3 0.000 0.000 0.000 0.000 0.000 0.208
## Compute the multiple effect sizes
( MT <- fit3@Model@def.function(.x.=x) )
##       MT1       MT2 
## 0.6024929 1.2049857
## Compute the jacobian for the 'defined parameters'
JAC <- lavaan:::lavJacobianD(func=fit3@Model@def.function, x=x)

## Compute the sampling covariance matrix using delta method
MT.VCOV <- JAC %*% VCOV %*% t(JAC)
## Add the variable names for ease of reference
dimnames(MT.VCOV) <- list(names(MT), names(MT))
MT.VCOV
##            MT1        MT2
## MT1 0.04040173 0.02234192
## MT2 0.02234192 0.04355176

Multiple-endpoint studies

###################################################
### code chunk number 7: ME1
###################################################
lower <- '11          
          5, 10'
## Convert a lower triangle data into a covariance matrix
Cov1 <- getCov(lower, diag=TRUE, names=c("x1", "x2"))

lower <- '12          
          6, 11'
## Convert a lower triangle data into a covariance matrix
Cov2 <- getCov(lower, diag=TRUE, names=c("x1", "x2"))

## Convert covariance matrices into a list
Cov <- list(Cov1, Cov2)

## Means for the two groups
Mean <- list(c(10,11), c(12,13))

## Sample sizes for the groups
N <- c(50, 50)

## Assuming homogeneity of covariance matrices by
## using the same labels: "sd1", "sd2", and "r"
model4 <- 'eta1 =~ c("sd1", "sd1")*x1
           eta2 =~ c("sd2", "sd2")*x2
           eta1 ~~ c("r", "r")*eta2
           x1 ~ c("m1_1", "m1_2")*1
           x2 ~ c("m2_1", "m2_2")*1     
           x1 ~~ 0*x1
           x2 ~~ 0*x2
           # Multiple endpoint effect size 1
           ME1 := (m1_2 - m1_1)/sd1
           # Multiple endpoint effect size 2
           ME2 := (m2_2 - m2_1)/sd2'

fit4 <- sem(model4, sample.cov=Cov, sample.mean=Mean, 
            sample.nobs=N, std.lv=TRUE, 
            sample.cov.rescale=FALSE)

## Obtain the free parameters in the model
( x <- fit4@Fit@x )
##  [1]  3.3911650  3.2403703  0.5005173 10.0000000 11.0000000  3.3911650
##  [7]  3.2403703  0.5005173 12.0000000 13.0000000
## Obtain the sampling covariance matrix of the parameter estimates
( VCOV <- vcov(fit4) )
##      sd1   sd2   r     m1_1  m2_1  sd1   sd2   r     m1_2  m2_2 
## sd1  0.057                                                      
## sd2  0.014 0.052                                                
## r    0.006 0.006 0.006                                          
## m1_1 0.000 0.000 0.000 0.230                                    
## m2_1 0.000 0.000 0.000 0.110 0.210                              
## sd1  0.057 0.014 0.006 0.000 0.000 0.057                        
## sd2  0.014 0.052 0.006 0.000 0.000 0.014 0.052                  
## r    0.006 0.006 0.006 0.000 0.000 0.006 0.006 0.006            
## m1_2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.230      
## m2_2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.110 0.210
## Compute the multivariate effect sizes
( ME <- fit4@Model@def.function(.x.=x) )
##       ME1       ME2 
## 0.5897678 0.6172134
## Compute the jacobian for 'defined parameters'
JAC <- lavaan:::lavJacobianD(func=fit4@Model@def.function, x=x)

## Compute the sampling covariance matrix using delta method
ME.VCOV <- JAC %*% VCOV %*% t(JAC)
## Add the variable names for ease of reference
dimnames(ME.VCOV) <- list(names(ME), names(ME))
ME.VCOV
##            ME1        ME2
## ME1 0.04173913 0.02047665
## ME2 0.02047665 0.04190476

Multiple treatment with multiple-endpoint studies

###################################################
### code chunk number 8: ME_MT
###################################################
## Covariance matrix of the control group
lower <- '11          
          5, 10'
## Convert a lower triangle data into a covariance matrix
Cov1 <- getCov(lower, diag=TRUE, names=c("x1", "x2"))

## Covariance matrix of the treatment group 1
lower <- '12          
          6, 11'
Cov2 <- getCov(lower, diag=TRUE, names=c("x1", "x2"))

## Covariance matrix of the treatment group 2
lower <- '13          
          7, 12'
Cov3 <- getCov(lower, diag=TRUE, names=c("x1", "x2"))

## Convert covariance matrices into a list
Cov <- list(Cov1, Cov2, Cov3)

## Means for the three groups
## 10 and 11 are the means for variables 1 and 2
Mean <- list(c(10,11), c(12,13), c(13,14))

## Sample sizes for the groups
N <- c(50, 50, 50)

## Assuming homogeneity of covariance matrices
model5 <- 'eta1 =~ c("sd1", "sd1", "sd1")*x1
           eta2 =~ c("sd2", "sd2", "sd2")*x2
           eta1 ~~ c("r", "r", "r")*eta2
           ## The subscripts 0, 1 and 2 represent the means
           ##  of the control and two  treatment groups
           x1 ~ c("m1_0", "m1_1", "m1_2")*1
           x2 ~ c("m2_0", "m2_1", "m2_2")*1
           ## The measurement errors are fixed at 0
           x1 ~~ 0*x1
           x2 ~~ 0*x2
           ## Multiple endpoint effect size 1 for treatment group 1
           ES1_1 := (m1_1 - m1_0)/sd1
           ## Multiple endpoint effect size 2 for treatment group 1
           ES2_1 := (m2_1 - m2_0)/sd2
           ## Multiple endpoint effect size 1 for treatment group 2
           ES1_2 := (m1_2 - m1_0)/sd1
           ## Multiple endpoint effect size 2 for treatment group 2
           ES2_2 := (m2_2 - m2_0)/sd2'

fit5 <- sem(model5, sample.cov=Cov, sample.mean=Mean, 
            sample.nobs=N, std.lv=TRUE, 
            sample.cov.rescale=FALSE)

## Obtain the free parameters in the model
( x <- fit5@Fit@x )
##  [1]  3.464102  3.316625  0.522233 10.000000 11.000000  3.464102  3.316625
##  [8]  0.522233 12.000000 13.000000  3.464102  3.316625  0.522233 13.000000
## [15] 14.000000
## Obtain the sampling covariance matrix of the parameter estimates
VCOV <- vcov(fit5)

## Compute the multivariate effect sizes
( ES <- fit5@Model@def.function(.x.=x) )
##     ES1_1     ES2_1     ES1_2     ES2_2 
## 0.5773503 0.6030227 0.8660254 0.9045340
## Compute the jacobian for 'defined parameters'
JAC <- lavaan:::lavJacobianD(func=fit5@Model@def.function, x=x)

## Compute the sampling covariance matrix using delta method
ES.VCOV <- JAC %*% VCOV %*% t(JAC)
## Add the variable names for ease of reference
dimnames(ES.VCOV) <- list(names(ES), names(ES))
ES.VCOV
##            ES1_1      ES2_1      ES1_2      ES2_2
## ES1_1 0.04111111 0.02120582 0.02166667 0.01091942
## ES2_1 0.02120582 0.04121212 0.01091942 0.02181818
## ES1_2 0.02166667 0.01091942 0.04250000 0.02160145
## ES2_2 0.01091942 0.02181818 0.02160145 0.04272727

Correlation matrix

###################################################
### code chunk number 9: CM1
###################################################
library("metaSEM")

## Sample correlation matrix
( C1 <- matrix(c(1,0.5,0.4,0.5,1,0.2,0.4,0.2,1), ncol=3,
               dimnames=list(c("x1","x2","x3"), 
                             c("x1","x2","x3"))) )
##     x1  x2  x3
## x1 1.0 0.5 0.4
## x2 0.5 1.0 0.2
## x3 0.4 0.2 1.0
###################################################
### code chunk number 10: CM2
###################################################
## Standard deviations
SD <- diag(c(1.2, 1.3, 1.4))

## Convert the correlation matrix to a covariance matrix
C2 <- SD %*% C1 %*% SD
dimnames(C2) <- list(c("x1","x2","x3"), 
                     c("x1","x2","x3"))
C2
##       x1    x2    x3
## x1 1.440 0.780 0.672
## x2 0.780 1.690 0.364
## x3 0.672 0.364 1.960
###################################################
### code chunk number 11: CM3
###################################################
## Calculate the sampling covariance matrix of
## the correlation matrix with n=50
asyCov(C2, n=50)
##             x2x1        x3x1        x3x2
## x2x1 0.011479608 0.001285708 0.005234689
## x3x1 0.001285708 0.014400010 0.007714304
## x3x2 0.005234689 0.007714304 0.018808189
###################################################
### code chunk number 12: CM4
###################################################
## Calculate the sampling covariance matrix of
## the covariance matrix with n=50
asyCov(C2, n=50, cor.analysis=FALSE)
##            x1x1        x2x1       x3x1        x2x2       x3x2        x3x3
## x1x1 0.08463759 0.045845312 0.03949821 0.024832537 0.02139507 0.018433045
## x2x1 0.04584531 0.062082103 0.02139445 0.053804376 0.02897177 0.009983792
## x3x1 0.03949821 0.021394450 0.06681716 0.011587980 0.03619270 0.053761187
## x2x2 0.02483254 0.053804376 0.01158798 0.116576318 0.02510828 0.005406977
## x3x2 0.02139507 0.028971765 0.03619270 0.025108276 0.07030471 0.029120349
## x3x3 0.01843305 0.009983792 0.05376119 0.005406977 0.02912035 0.156801459

Chapter 4: Univariate Meta-Analysis

Odds ratio of atrial fibril lation between bisphosphonate and non-bisphosphonate users

###################################################
### code chunk number 2: Mak1
###################################################
## Display the full dataset
Mak09
##                Study type AF.BP Tot.BP AF.non.BP Tot.non.BP          yi
## 1       Black (2007)  RCT    94   3862        73       3852  0.25575041
## 2    Cummings (2007)  RCT    81   3236        71       3223  0.13081795
## 3       Karam (2007)  RCT   189  10018        94       5048  0.01331036
## 4       Lyles (2007)  RCT    29   1054        27       1057  0.07632515
## 5 Papapoulous (2008)  RCT    57   6830        18       1924 -0.11525781
## 6  Abrahamsen (2009)  Obs   797  14302      1280      28731  0.23558193
## 7    Heckbert (2008)  Obs    47     87       672       1598  0.48188404
## 8    Sorensen (2008)  Obs   724   3862     12862      77643  0.15018558
##            vi age.mean study.duration
## 1 0.024866941     73.0            3.0
## 2 0.027064402     69.3            3.6
## 3 0.016232900     73.5             NA
## 4 0.073466279     74.5            5.0
## 5 0.073771719     66.9            2.5
## 6 0.002146430     74.3           10.0
## 7 0.048844605     72.7            3.0
## 8 0.001793075     76.1            6.0
###################################################
### code chunk number 3: Mak2 
###################################################
mak1 <- meta(y=yi, v=vi, data=Mak09)
summary(mak1)
## 
## Call:
## meta(y = yi, v = vi, data = Mak09)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value
## Intercept1  1.8082e-01  2.8942e-02  1.2410e-01  2.3755e-01  6.2477
## Tau2_1_1    1.6561e-10  2.2166e-03 -4.3444e-03  4.3444e-03  0.0000
##             Pr(>|z|)    
## Intercept1 4.165e-10 ***
## Tau2_1_1           1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 7.159836
## Degrees of freedom of the Q statistic: 7
## P value of the Q statistic: 0.4124299
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)        0
## 
## Number of studies (or clusters): 8
## Number of observed statistics: 8
## Number of estimated parameters: 2
## Degrees of freedom: 6
## -2 log likelihood: -10.26621 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
#### Alternative specification without storing the results
summary( meta(y=yi, v=vi, data=Mak09) )
## 
## Call:
## meta(y = yi, v = vi, data = Mak09)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value
## Intercept1  1.8082e-01  2.8942e-02  1.2410e-01  2.3755e-01  6.2477
## Tau2_1_1    1.6561e-10  2.2166e-03 -4.3444e-03  4.3444e-03  0.0000
##             Pr(>|z|)    
## Intercept1 4.165e-10 ***
## Tau2_1_1           1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 7.159836
## Degrees of freedom of the Q statistic: 7
## P value of the Q statistic: 0.4124299
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)        0
## 
## Number of studies (or clusters): 8
## Number of observed statistics: 8
## Number of estimated parameters: 2
## Degrees of freedom: 6
## -2 log likelihood: -10.26621 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 5: Mak4
###################################################
( mak2 <- summary(meta(y=yi, v=vi, data=Mak09, RE.constraints=0)) )
## 
## Call:
## meta(y = yi, v = vi, data = Mak09, RE.constraints = 0)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##            Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## Intercept1 0.180822  0.028748 0.124477 0.237167  6.2899 3.177e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 7.159836
## Degrees of freedom of the Q statistic: 7
## P value of the Q statistic: 0.4124299
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)        0
## 
## Number of studies (or clusters): 8
## Number of observed statistics: 8
## Number of estimated parameters: 1
## Degrees of freedom: 7
## -2 log likelihood: -10.26621 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
#### Alternative specification
( mak2 <- summary(meta(y=yi, v=vi, data=Mak09,
                   RE.constraints=matrix(0, ncol=1, nrow=1))) )
## 
## Call:
## meta(y = yi, v = vi, data = Mak09, RE.constraints = matrix(0, 
##     ncol = 1, nrow = 1))
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##            Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## Intercept1 0.180822  0.028748 0.124477 0.237167  6.2899 3.177e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 7.159836
## Degrees of freedom of the Q statistic: 7
## P value of the Q statistic: 0.4124299
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)        0
## 
## Number of studies (or clusters): 8
## Number of observed statistics: 8
## Number of estimated parameters: 1
## Degrees of freedom: 7
## -2 log likelihood: -10.26621 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 7: Mak6
###################################################
## Get the estimate and its 95% CI
( Est <- mak2$coefficients["Intercept1", 
                          c("Estimate","lbound","ubound")] )
##             Estimate    lbound    ubound
## Intercept1 0.1808222 0.1244772 0.2371673
## Convert them into odds ratio
exp(Est)
##            Estimate   lbound   ubound
## Intercept1 1.198202 1.132556 1.267653

Correlation between organizational commitment and salesperson job performance

###################################################
### code chunk number 8: Jaramillo1
###################################################
head(Jaramillo05)
##                        Author Sample_size    Sales Country IDV
## 1         Aryee et al. (2002)         179    mixed   India  48
## 2 Balfour and Wechsler (1991)         232 nonsales     USA  91
## 3     Bashaw and Grant (1994)         560    sales     USA  91
## 4             Benkhoff (1997)         181    sales Germany  67
## 5         Brett et al. (1995)         156    sales     USA  91
## 6         Brett et al. (1995)         180    sales     USA  91
##           OC_scale OC_alpha JP_alpha    r         r_v
## 1 Porter or Mowday     0.87     0.89 0.02 0.005582124
## 2            other     0.82       NA 0.12 0.004187101
## 3 Porter or Mowday     0.83     0.76 0.09 0.001756903
## 4 Porter or Mowday       NA     1.00 0.20 0.005091713
## 5 Porter or Mowday     0.83       NA 0.08 0.006328468
## 6 Porter or Mowday     0.83       NA 0.04 0.005537792
###################################################
### code chunk number 9: Jaramillo2 (eval = FALSE)
###################################################
## z <- with( Jaramillo05, 0.5*log((1+r)/(1-r)) )
## z.v <- with( Jaramillo05, 1/(Sample.size-3) )


###################################################
### code chunk number 10: Jaramillo3
###################################################
summary( meta(y=r, v=r_v, data=Jaramillo05) )
## 
## Call:
## meta(y = r, v = r_v, data = Jaramillo05)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##             Estimate Std.Error    lbound    ubound z value  Pr(>|z|)    
## Intercept1 0.1866221 0.0193303 0.1487354 0.2245088  9.6544 < 2.2e-16 ***
## Tau2_1_1   0.0170341 0.0041352 0.0089292 0.0251389  4.1193 3.801e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 339.3886
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.8144
## 
## Number of studies (or clusters): 61
## Number of observed statistics: 61
## Number of estimated parameters: 2
## Degrees of freedom: 59
## -2 log likelihood: -55.44225 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 12: Jaramillo5
###################################################
summary( meta(y=r, v=r_v, data=Jaramillo05, intervals.type="LB") )
## 
## Call:
## meta(y = r, v = r_v, data = Jaramillo05, intervals.type = "LB")
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##            Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## Intercept1 0.186622        NA 0.147975 0.225074      NA       NA
## Tau2_1_1   0.017034        NA 0.010596 0.026388      NA       NA
## 
## Q statistic on the homogeneity of effect sizes: 339.3886
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Heterogeneity indices (I2) and their 95% likelihood-based CIs:
##                               lbound Estimate ubound
## Intercept1: I2 (Q statistic) 0.73193  0.81438 0.8404
## 
## Number of studies (or clusters): 61
## Number of observed statistics: 61
## Number of estimated parameters: 2
## Degrees of freedom: 59
## -2 log likelihood: -55.44225 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 14: Jaramillo7 
###################################################
## Center IND: scale(IND, scale=FALSE)
##  scale=TRUE: standardize the variable
##  scale=FALSE: not standardize the variable
model0 <- meta(y=r, v=r_v, x=scale(IDV, scale=FALSE), 
              data=Jaramillo05)

summary(model0)
## 
## Call:
## meta(y = r, v = r_v, x = scale(IDV, scale = FALSE), data = Jaramillo05)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value
## Intercept1  0.18596328  0.01903036  0.14866446  0.22326210  9.7719
## Slope1_1   -0.00132138  0.00097695 -0.00323617  0.00059341 -1.3526
## Tau2_1_1    0.01634054  0.00402026  0.00846097  0.02422010  4.0645
##             Pr(>|z|)    
## Intercept1 < 2.2e-16 ***
## Slope1_1      0.1762    
## Tau2_1_1   4.813e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 339.3886
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Explained variances (R2):
##                            y1
## Tau2 (no predictor)    0.0170
## Tau2 (with predictors) 0.0163
## R2                     0.0407
## 
## Number of studies (or clusters): 61
## Number of observed statistics: 61
## Number of estimated parameters: 3
## Degrees of freedom: 58
## -2 log likelihood: -57.24231 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 16: Jaramillo9 
###################################################
model1 <- meta(y=r, v=r_v, x=cbind(OC_alpha, JP_alpha), 
               data=Jaramillo05,
               model.name="Unequal coefficients")
summary(model1)
## 
## Call:
## meta(y = r, v = r_v, x = cbind(OC_alpha, JP_alpha), data = Jaramillo05, 
##     model.name = "Unequal coefficients")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
## Intercept1 -0.5755371  0.5015222 -1.5585026  0.4074284 -1.1476 0.2511418
## Slope1_1    0.1311033  0.4587248 -0.7679807  1.0301873  0.2858 0.7750317
## Slope1_2    0.8044184  0.4303887 -0.0391280  1.6479648  1.8691 0.0616158
## Tau2_1_1    0.0187312  0.0056537  0.0076503  0.0298122  3.3131 0.0009226
##               
## Intercept1    
## Slope1_1      
## Slope1_2   .  
## Tau2_1_1   ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 257.7036
## Degrees of freedom of the Q statistic: 34
## P value of the Q statistic: 0
## 
## Explained variances (R2):
##                            y1
## Tau2 (no predictor)    0.0170
## Tau2 (with predictors) 0.0187
## R2                     0.0000
## 
## Number of studies (or clusters): 35
## Number of observed statistics: 35
## Number of estimated parameters: 4
## Degrees of freedom: 31
## -2 log likelihood: -31.01205 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 18: Jaramillo11a
###################################################
( constraint <- matrix(c("0*Slope_equal", "0*Slope_equal"), 
                       nrow=1, ncol=2) )
##      [,1]            [,2]           
## [1,] "0*Slope_equal" "0*Slope_equal"
###################################################
### code chunk number 19: Jaramillo11b
###################################################
model2 <- meta(y=r, v=r_v, x=cbind(OC_alpha, JP_alpha), 
               data=Jaramillo05, coef.constraints=constraint, 
               model.name="Equal coefficients")
summary(model2)
## 
## Call:
## meta(y = r, v = r_v, x = cbind(OC_alpha, JP_alpha), data = Jaramillo05, 
##     coef.constraints = constraint, model.name = "Equal coefficients")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
## Intercept1  -0.6036671  0.5074490 -1.5982488  0.3909146 -1.1896 0.2341991
## Slope_equal  0.4863010  0.2953055 -0.0924871  1.0650892  1.6468 0.0996048
## Tau2_1_1     0.0193743  0.0058239  0.0079596  0.0307890  3.3267 0.0008789
##                
## Intercept1     
## Slope_equal .  
## Tau2_1_1    ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 257.7036
## Degrees of freedom of the Q statistic: 34
## P value of the Q statistic: 0
## 
## Explained variances (R2):
##                            y1
## Tau2 (no predictor)    0.0170
## Tau2 (with predictors) 0.0194
## R2                     0.0000
## 
## Number of studies (or clusters): 35
## Number of observed statistics: 35
## Number of estimated parameters: 3
## Degrees of freedom: 32
## -2 log likelihood: -30.01779 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 21: Jaramillo11
###################################################
anova(model1, model2)
##                   base         comparison ep  minus2LL df       AIC
## 1 Unequal coefficients               <NA>  4 -31.01205 31 -93.01205
## 2 Unequal coefficients Equal coefficients  3 -30.01779 32 -94.01779
##      diffLL diffdf        p
## 1        NA     NA       NA
## 2 0.9942617      1 0.318703
###################################################
### code chunk number 23: Jaramillo14
###################################################
table(Jaramillo05$Sales)
## 
##    mixed nonsales    sales 
##        6       27       28
sales <- ifelse(Jaramillo05$Sales=="sales", yes=1, no=0)
nonsales <- ifelse(Jaramillo05$Sales=="nonsales", yes=1, no=0)
mixed <- ifelse(Jaramillo05$Sales=="mixed", yes=1, no=0)


###################################################
### code chunk number 24: Jaramillo15a
###################################################
( startvalues <- matrix(c("0*Slope1_1", "0*Slope1_2", 
                          "0*Slope1_3"), nrow=1, ncol=3) )
##      [,1]         [,2]         [,3]        
## [1,] "0*Slope1_1" "0*Slope1_2" "0*Slope1_3"
###################################################
### code chunk number 25: Jaramillo15b
###################################################
model3 <- meta(y=r, v=r_v, x=cbind(sales, mixed, nonsales), 
               data=Jaramillo05, coef.constraints=startvalues,
               intercept.constraints=matrix(0, ncol=1, nrow=1),               
               model.name="Indicator variables")
summary(model3)
## 
## Call:
## meta(y = r, v = r_v, x = cbind(sales, mixed, nonsales), data = Jaramillo05, 
##     intercept.constraints = matrix(0, ncol = 1, nrow = 1), coef.constraints = startvalues, 
##     model.name = "Indicator variables")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##           Estimate Std.Error    lbound    ubound z value  Pr(>|z|)    
## Slope1_1 0.2282976 0.0275950 0.1742125 0.2823828  8.2732 2.220e-16 ***
## Slope1_2 0.1465871 0.0632760 0.0225685 0.2706057  2.3166   0.02052 *  
## Slope1_3 0.1519563 0.0279388 0.0971973 0.2067152  5.4389 5.361e-08 ***
## Tau2_1_1 0.0157300 0.0038524 0.0081793 0.0232806  4.0831 4.443e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 339.3886
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Explained variances (R2):
##                            y1
## Tau2 (no predictor)    0.0170
## Tau2 (with predictors) 0.0157
## R2                     0.0766
## 
## Number of studies (or clusters): 61
## Number of observed statistics: 61
## Number of estimated parameters: 4
## Degrees of freedom: 57
## -2 log likelihood: -59.55628 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 27: Jaramillo17
###################################################
model4 <- meta(y=r, v=r_v, data=Jaramillo05)
anova(model3, model4)
##                  base            comparison ep  minus2LL df       AIC
## 1 Indicator variables                  <NA>  4 -59.55628 57 -173.5563
## 2 Indicator variables Meta analysis with ML  2 -55.44225 59 -173.4423
##     diffLL diffdf         p
## 1       NA     NA        NA
## 2 4.114024      2 0.1278354

Chapter 5: Multivariate Meta-Analysis

BCG vaccine for preventing tuberculosis

###################################################
### code chunk number 2: BCG1
###################################################
#### Load the metaSEM library
library("metaSEM")
## Display the dataset
head(BCG)
##   Trial               Author Year  VD   VWD NVD  NVWD Latitude Allocation
## 1     1              Aronson 1948   4   119  11   128       44     random
## 2     2     Ferguson & Simes 1949   6   300  29   274       55     random
## 3     3      Rosenthal et al 1960   3   228  11   209       42     random
## 4     4    Hart & Sutherland 1977  62 13536 248 12619       52     random
## 5     5 Frimodt-Moller et al 1973  33  5036  47  5761       13  alternate
## 6     6      Stein & Aronson 1953 180  1361 372  1079       44  alternate
##        ln_OR     v_ln_OR  ln_Odd_V ln_Odd_NV  v_ln_Odd_V cov_V_NV
## 1 -0.9386941 0.357124952 -3.392829 -2.454135 0.258403361        0
## 2 -1.6661907 0.208132394 -3.912023 -2.245832 0.170000000        0
## 3 -1.3862944 0.433413078 -4.330733 -2.944439 0.337719298        0
## 4 -1.4564435 0.020314413 -5.385974 -3.929530 0.016202909        0
## 5 -0.2191411 0.051951777 -5.027860 -4.808719 0.030501601        0
## 6 -0.9581220 0.009905266 -2.023018 -1.064896 0.006290309        0
##   v_ln_Odd_NV
## 1 0.098721591
## 2 0.038132394
## 3 0.095693780
## 4 0.004111504
## 5 0.021450177
## 6 0.003614956
###################################################
### code chunk number 3: BCG2
###################################################
## Covariance between the effect size is 0.
bcg1 <- meta(y=cbind(ln_Odd_V, ln_Odd_NV), 
             v=cbind(v_ln_Odd_V, cov_V_NV, v_ln_Odd_NV),
             data=BCG, model.name="Random effects model") 
summary(bcg1)
## 
## Call:
## meta(y = cbind(ln_Odd_V, ln_Odd_NV), v = cbind(v_ln_Odd_V, cov_V_NV, 
##     v_ln_Odd_NV), data = BCG, model.name = "Random effects model")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##            Estimate Std.Error   lbound   ubound  z value Pr(>|z|)    
## Intercept1 -4.83374   0.34020 -5.50052 -4.16697 -14.2086  < 2e-16 ***
## Intercept2 -4.09597   0.43474 -4.94806 -3.24389  -9.4216  < 2e-16 ***
## Tau2_1_1    1.43137   0.58304  0.28863  2.57411   2.4550  0.01409 *  
## Tau2_2_1    1.75733   0.72425  0.33782  3.17683   2.4264  0.01525 *  
## Tau2_2_2    2.40733   0.96742  0.51123  4.30344   2.4884  0.01283 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 5270.386
## Degrees of freedom of the Q statistic: 24
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.9887
## Intercept2: I2 (Q statistic)   0.9955
## 
## Number of studies (or clusters): 13
## Number of observed statistics: 26
## Number of estimated parameters: 5
## Degrees of freedom: 21
## -2 log likelihood: 66.17587 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 5: BCG4
###################################################
bcg2 <- meta(y=cbind(ln_Odd_V, ln_Odd_NV), data=BCG,
             v=cbind(v_ln_Odd_V, cov_V_NV, v_ln_Odd_NV),
             intercept.constraints=c("0*Intercept","0*Intercept"),
             model.name="Equal intercepts") 
summary(bcg2)
## 
## Call:
## meta(y = cbind(ln_Odd_V, ln_Odd_NV), v = cbind(v_ln_Odd_V, cov_V_NV, 
##     v_ln_Odd_NV), data = BCG, intercept.constraints = c("0*Intercept", 
##     "0*Intercept"), model.name = "Equal intercepts")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##            Estimate Std.Error    lbound    ubound  z value Pr(>|z|)    
## Intercept -5.375013  0.458434 -6.273527 -4.476498 -11.7247  < 2e-16 ***
## Tau2_1_1   1.700312  0.855140  0.024269  3.376355   1.9883  0.04677 *  
## Tau2_2_1   2.455619  1.311343 -0.114566  5.025804   1.8726  0.06112 .  
## Tau2_2_2   4.108667  1.990284  0.207783  8.009551   2.0644  0.03898 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 5270.386
## Degrees of freedom of the Q statistic: 24
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.9905
## Intercept2: I2 (Q statistic)   0.9974
## 
## Number of studies (or clusters): 13
## Number of observed statistics: 26
## Number of estimated parameters: 4
## Degrees of freedom: 22
## -2 log likelihood: 77.05836 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 7: BCG6
###################################################
anova(bcg1, bcg2)
##                   base       comparison ep minus2LL df      AIC  diffLL
## 1 Random effects model             <NA>  5 66.17587 21 24.17587      NA
## 2 Random effects model Equal intercepts  4 77.05836 22 33.05836 10.8825
##   diffdf            p
## 1     NA           NA
## 2      1 0.0009707732
###################################################
### code chunk number 9: BCG7
###################################################
## Extract the coefficient table from the summary
Est <- summary(bcg1)$coefficients

## Only select the first 2 rows and the columns
## related to estimate, lbound, and ubound
## Convert them into odds
exp( Est[1:2, c("Estimate", "lbound", "ubound") ] )
##               Estimate      lbound     ubound
## Intercept1 0.007956678 0.004084650 0.01549918
## Intercept2 0.016639521 0.007097173 0.03901183
###################################################
### code chunk number 10: BCG8
###################################################
## Extract the fixed effects 
( fixed <- coef(bcg1, select="fixed") )
## Intercept1 Intercept2 
##  -4.833744  -4.095975
## Extract the sampling covariance matrix on the estimates
( omega <- vcov(bcg1)[c("Intercept1","Intercept2"), 
                      c("Intercept1","Intercept2")] )
##            Intercept1 Intercept2
## Intercept1  0.1157346  0.1359696
## Intercept2  0.1359696  0.1890030
## Calculate the logarithm on the odds ratio
( log_OR <- fixed[1] - fixed[2] )
## Intercept1 
## -0.7377691
## Calculate the standard error on log_OR
( se_log_OR <- sqrt(omega[1,1]+omega[2,2]-2*omega[2,1]) )
## [1] 0.1811031
###################################################
### code chunk number 11: BCG9 
###################################################
summary( meta(y=ln_OR, v=v_ln_OR, data=BCG) )
## 
## Call:
## meta(y = ln_OR, v = v_ln_OR, data = BCG)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
## Intercept1 -0.7419668  0.1786345 -1.0920841 -0.3918496 -4.1535 3.274e-05
## Tau2_1_1    0.3024566  0.1566294 -0.0045315  0.6094446  1.9310   0.05348
##               
## Intercept1 ***
## Tau2_1_1   .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 163.1649
## Degrees of freedom of the Q statistic: 12
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.9123
## 
## Number of studies (or clusters): 13
## Number of observed statistics: 13
## Number of estimated parameters: 2
## Degrees of freedom: 11
## -2 log likelihood: 26.14552 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 13: BCG11
###################################################
bcg3 <- meta(y=cbind(ln_Odd_V, ln_Odd_NV), data=BCG,
             v=cbind(v_ln_Odd_V, cov_V_NV, v_ln_Odd_NV),
             RE.constraints=matrix(c("0.1*Tau2_Eq","0*Tau2_2_1",
                                     "0*Tau2_2_1","0.1*Tau2_Eq"),
                                   ncol=2, nrow=2),
             model.name="Equal variances")
summary(bcg3)
## 
## Call:
## meta(y = cbind(ln_Odd_V, ln_Odd_NV), v = cbind(v_ln_Odd_V, cov_V_NV, 
##     v_ln_Odd_NV), data = BCG, RE.constraints = matrix(c("0.1*Tau2_Eq", 
##     "0*Tau2_2_1", "0*Tau2_2_1", "0.1*Tau2_Eq"), ncol = 2, nrow = 2), 
##     model.name = "Equal variances")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##            Estimate Std.Error   lbound   ubound  z value  Pr(>|z|)    
## Intercept1 -4.83659   0.39241 -5.60569 -4.06748 -12.3254 < 2.2e-16 ***
## Intercept2 -4.08256   0.38911 -4.84520 -3.31991 -10.4920 < 2.2e-16 ***
## Tau2_Eq     1.91972   0.73759  0.47407  3.36536   2.6027  0.009249 ** 
## Tau2_2_1    1.76544   0.73886  0.31730  3.21358   2.3894  0.016875 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 5270.386
## Degrees of freedom of the Q statistic: 24
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.9915
## Intercept2: I2 (Q statistic)   0.9944
## 
## Number of studies (or clusters): 13
## Number of observed statistics: 26
## Number of estimated parameters: 4
## Degrees of freedom: 22
## -2 log likelihood: 71.24705 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 15: BCG13
###################################################
anova(bcg1, bcg3)
##                   base      comparison ep minus2LL df      AIC   diffLL
## 1 Random effects model            <NA>  5 66.17587 21 24.17587       NA
## 2 Random effects model Equal variances  4 71.24705 22 27.24705 5.071183
##   diffdf          p
## 1     NA         NA
## 2      1 0.02432678
###################################################
### code chunk number 17: BCG15
###################################################
( T2 <- vec2symMat(coef(bcg1, select="random")) )
##          [,1]     [,2]
## [1,] 1.431371 1.757327
## [2,] 1.757327 2.407333
###################################################
### code chunk number 18: BCG16
###################################################
cov2cor(T2)
##           [,1]      [,2]
## [1,] 1.0000000 0.9466912
## [2,] 0.9466912 1.0000000
###################################################
### code chunk number 19: BCG17
###################################################
bcg.cor <- cov2cor(T2)[2,1]

###################################################
### code chunk number 20: BCG18 
###################################################
plot(bcg1, xlim=c(-8,0), ylim=c(-8,0))

###################################################
### code chunk number 21: BCG19
###################################################
## ## Load the metafor package
library("metafor")
plot(bcg1, xlim=c(-8,0), ylim=c(-8,0), diag.panel=TRUE)
## Forest plot for the vaccinated group
forest( rma(yi=ln_Odd_V, vi=v_ln_Odd_V, method="ML", data=BCG) )
title("Forest plot for the vaccinated group")
## Forest plot for the non-vaccinated group
forest( rma(yi=ln_Odd_NV, vi=v_ln_Odd_NV, method="ML", data=BCG) )
title("Forest plot for the non-vaccinated group")

Standardized mean differences between males and females on life satisfaction and life control

###################################################
### code chunk number 22: wvs1
###################################################
## Display the dataset
head(wvs94a)
##     country      lifesat    lifecon lifesat_var   inter_cov lifecon_var
## 1 Argentina -0.032092632 0.05760759 0.004043047 0.001407526 0.004158255
## 2   Austria  0.080096043 0.00889294 0.002887624 0.000933741 0.002893639
## 3   Belarus  0.041978879 0.07408653 0.004009652 0.001335892 0.004010990
## 4   Belgium  0.007754997 0.12799510 0.001456405 0.000405010 0.001513248
## 5    Brazil  0.148138118 0.18210572 0.002266132 0.000789084 0.002290195
## 6   Britain  0.020047940 0.04445499 0.002723780 0.001185755 0.002746361
##     gnp
## 1  2370
## 2  4900
## 3  3110
## 4 15540
## 5  2680
## 6 16100
###################################################
### code chunk number 23: wvs2 
###################################################
## Random-effects model
wvs1 <- meta(y=cbind(lifesat, lifecon),
             v=cbind(lifesat_var, inter_cov, lifecon_var), 
             data=wvs94a, model.name="Random effects model")
summary(wvs1)
## 
## Call:
## meta(y = cbind(lifesat, lifecon), v = cbind(lifesat_var, inter_cov, 
##     lifecon_var), data = wvs94a, model.name = "Random effects model")
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value
## Intercept1  0.00134985  0.01385628 -0.02580796  0.02850767  0.0974
## Intercept2  0.06882575  0.01681962  0.03585990  0.10179160  4.0920
## Tau2_1_1    0.00472726  0.00176156  0.00127465  0.00817986  2.6836
## Tau2_2_1    0.00393437  0.00168706  0.00062779  0.00724094  2.3321
## Tau2_2_2    0.00841361  0.00253727  0.00344064  0.01338657  3.3160
##             Pr(>|z|)    
## Intercept1 0.9223943    
## Intercept2 4.277e-05 ***
## Tau2_1_1   0.0072844 ** 
## Tau2_2_1   0.0196962 *  
## Tau2_2_2   0.0009131 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 250.0303
## Degrees of freedom of the Q statistic: 82
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.6129
## Intercept2: I2 (Q statistic)   0.7345
## 
## Number of studies (or clusters): 42
## Number of observed statistics: 84
## Number of estimated parameters: 5
## Degrees of freedom: 79
## -2 log likelihood: -161.9216 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
###################################################
### code chunk number 25: wvs4
###################################################
## Extract the variance component of the random effects
( T2 <- vec2symMat(coef(wvs1, select="random")) )
##             [,1]        [,2]
## [1,] 0.004727256 0.003934366
## [2,] 0.003934366 0.008413607
###################################################
### code chunk number 26: wvs5
###################################################
## Convert the covariance matrix to a correlation matrix
cov2cor(T2)
##           [,1]      [,2]
## [1,] 1.0000000 0.6238485
## [2,] 0.6238485 1.0000000
###################################################
### code chunk number 27: wvs6
###################################################
wvs.cor <- cov2cor(T2)[2,1]


###################################################
### code chunk number 28: wvs7
###################################################
par(mfrow=c(1,1))
plot(wvs1, axis.labels=c("SMD on life satisfaction", 
                        "SMD on life control"),
    study.ellipse.plot=FALSE,
    xlim=c(-0.3, 0.2), ylim=c(-0.3,0.4))