Simulation meta analysis

Simulation for meta analysis

Goal: simulate data of two indepentent BLUP predictors of a random quantity and evaluate alternative meta-analysis weight combinations.

Simulation

Start by setting up inputs: variance of random variable and variances of estimators (have to be smaller than the variance of the variable). Number of simulations. Compute PEV

vg<-2
v1<-0.5
v2<-0.3
nsim=1000000

pev1<-vg-v1
pev2<-vg-v2
pev1

[1] 1.5

pev2

[1] 1.7

Build V-CV matrix of the three quantities (variable and estimates). Compute conditional expectations of estimates conditional on the variable value (muc). Compute conditional variance too (vc).

Simulate estimates (gh) using conditional variance and expectations.

Compute empirical variances of simulated variable and estimates (EV). Check that values are consistent with inputs (SS). Also check empirical PEVs.

g<-rnorm(nsim,sd=sqrt(vg))
SS<-rbind(c(vg,v1,v2),c(v1,v1,0),c(v2,0,v2))
SS

     [,1] [,2] [,3]
[1,]  2.0  0.5  0.3
[2,]  0.5  0.5  0.0
[3,]  0.3  0.0  0.3

muc<-  SS[2:3,1]/SS[1,1]
muc

[1] 0.25 0.15

mn<-g%*%matrix(muc,1,2)

matrix(muc,2,1)

     [,1]
[1,] 0.25
[2,] 0.15

vc<-SS[2:3,2:3]-SS[2:3,1,drop=F]%*%SS[1,2:3,drop=F]/SS[1,1]


bs<-matrix(rnorm(2*nsim),nsim,2)
gh<-mn+bs%*%chol(vc)
dim(gh)

[1] 1000000       2

cov(gh)

              [,1]          [,2]
[1,]  0.4994912296 -0.0001359151
[2,] -0.0001359151  0.3005480512

EV<-var(cbind(g,gh))
EV

          g                            
g 1.9994331  0.4987812167  0.3000270946
  0.4987812  0.4994912296 -0.0001359151
  0.3000271 -0.0001359151  0.3005480512

SS

     [,1] [,2] [,3]
[1,]  2.0  0.5  0.3
[2,]  0.5  0.5  0.0
[3,]  0.3  0.0  0.3

EV[1,1]-EV[2,2]

[1] 1.499942

EV[1,1]-EV[3,3]

[1] 1.698885

pev1

[1] 1.5

pev2

[1] 1.7

Meta analysis

Set up weights, in this case, 1/pev, test other values.

Compute meta-analysis predictor (gm).

Compare PEV, accuracy and slope of estimates compared to random variable

w1=1/pev1
w2=1/pev2


gm<-(w1*gh[,1]+w2*gh[,2])/(w1+w2)


#exmpirical PEV
var(gm-g)

[1] 1.395143

var(gh[,1]-g)

[1] 1.501362

var(gh[,2]-g)

[1] 1.699927

#theoretical PEV
pev1

[1] 1.5

pev2

[1] 1.7

pevm<-vg - (w1^2*v1+w2^2*v2)/((w1+w2)^2)
pevm #WRONG!

[1] 1.792969

var(gm)

[1] 0.2069404

(w1^2*v1+w2^2*v2)/((w1+w2)^2) #Right

[1] 0.2070312

cov(gm,g) #differemt from var(gm)!!!

[1] 0.4056152

pevm2<-vg + (((w1^2)*v1+(w2^2)*v2)/((w1+w2)^2)) - (2*(w1*v1+w2*v2)/(w1+w2))
pevm2 #RIGHT!!!

[1] 1.394531

#empirical accuracies
cor(gm,g)

[1] 0.6305771

cor(g,gh[,1])

[1] 0.4991059

cor(g,gh[,1])

[1] 0.4991059

coefficients(lm(g~gm)) #out of scale!

 (Intercept)           gm 
0.0009623558 1.9600582584

coefficients(lm(g~gh[,1]))

(Intercept)     gh[, 1] 
0.001292536 0.998578528

coefficients(lm(g~gh[,1]))

(Intercept)     gh[, 1] 
0.001292536 0.998578528

#theoretical accuracy from PEV
sqrt(1-pevm2/vg) #it does not check!

[1] 0.550213

Conclusion:

We need to derive the optimal weights using these equations. (maximize correlation and minimize PEVm)
We need to work on the scale of MA estimates