metahdep.R
In metahdep: Hierarchical Dependence in Meta-Analysis

### R code from vignette source 'metahdep.Rnw'

###################################################
### code chunk number 1: metahdep.Rnw:155-164 (eval = FALSE)
###################################################
## library(AnnotationDbi)
## library(hgu133a.db); library(hgu133b.db); library(hgu95a.db); library(hgu95av2.db)
## EID.hgu133a <- unlist(mget(mappedkeys(hgu133aENTREZID),hgu133aENTREZID))
## EID.hgu133b <- unlist(mget(mappedkeys(hgu133bENTREZID),hgu133bENTREZID))
## EID.hgu95a <- unlist(mget(mappedkeys(hgu95aENTREZID),hgu95aENTREZID))
## EID.hgu95av2 <- unlist(mget(mappedkeys(hgu95av2ENTREZID),hgu95av2ENTREZID))
## chip <- c(rep("hgu133a", length(EID.hgu133a)), rep("hgu133b", length(EID.hgu133b)), rep("hgu95a", length(EID.hgu95a)), rep("hgu95av2", length(EID.hgu95av2)))
## EID <- c(EID.hgu133a, EID.hgu133b, EID.hgu95a, EID.hgu95av2)
## newnames <- data.frame(chip=chip, old.name=names(EID), new.name=EID)


###################################################
### code chunk number 2: metahdep.Rnw:174-181 (eval = FALSE)
###################################################
## set.seed(1234)
## use.new.name <- sample(unique(newnames$new.name),1000)
## use.hgu133a.gn <- names(EID.hgu133a[is.element(EID.hgu133a,use.new.name)])
## use.hgu133b.gn <- names(EID.hgu133b[is.element(EID.hgu133b,use.new.name)])
## use.hgu95a.gn <- names(EID.hgu95a[is.element(EID.hgu95a,use.new.name)])
## use.hgu95av2.gn <- names(EID.hgu95av2[is.element(EID.hgu95av2,use.new.name)])
## HGU.newnames <- newnames[is.element(EID,use.new.name),]


###################################################
### code chunk number 3: metahdep.Rnw:242-244 (eval = FALSE)
###################################################
## library(ALLMLL); data(MLL.A)
## ES.exp1 <- getPLM.es(abatch=MLL.A, trt1=list(c(1:5),c(6:10)), trt2=list(c(11:15),c(16:20)), sub.gn=use.hgu133a.gn, covariates=data.frame(Tissue=c(0,1)), dep.grp=1 )


###################################################
### code chunk number 4: metahdep.Rnw:259-262 (eval = FALSE)
###################################################
## library(SpikeInSubset)
## data(spikein95)
## ES.exp2 <- getPLM.es(abatch=spikein95, trt1=1:3, trt2=4:6, sub.gn=use.hgu95a.gn, covariates=data.frame(Tissue=0), dep.grp=1 )


###################################################
### code chunk number 5: metahdep.Rnw:271-273 (eval = FALSE)
###################################################
## data(Dilution)
## ES.exp3 <- getPLM.es(abatch=Dilution, trt1=1:2, trt2=3:4, sub.gn=use.hgu95av2.gn, covariates=data.frame(Tissue=0), dep.grp=2 )


###################################################
### code chunk number 6: metahdep.Rnw:282-284 (eval = FALSE)
###################################################
## data(MLL.B)
## ES.exp4 <- getPLM.es(abatch=MLL.B, trt1=1:9, trt2=10:20, sub.gn=use.hgu133b.gn, covariates=data.frame(Tissue=1), dep.grp=3 )


###################################################
### code chunk number 7: metahdep.Rnw:294-295 (eval = FALSE)
###################################################
## HGU.DifExp.list <- list(ES.exp1, ES.exp2, ES.exp3, ES.exp4)


###################################################
### code chunk number 8: metahdep.Rnw:313-314 (eval = FALSE)
###################################################
## HGU.prep.list <- metahdep.format(HGU.DifExp.list, HGU.newnames)


###################################################
### code chunk number 9: metahdep.Rnw:348-353
###################################################
library(metahdep)
data(HGU.newnames)
data(HGU.prep.list)
set.seed(123)
gene.subset <- sample(HGU.newnames$new.name, 50)


###################################################
### code chunk number 10: metahdep.Rnw:358-359
###################################################
FEMA <- metahdep(HGU.prep.list, method="FEMA", genelist=gene.subset, center.X=TRUE)


###################################################
### code chunk number 11: metahdep.Rnw:366-368
###################################################
REMA <- metahdep(HGU.prep.list, method="REMA", genelist=gene.subset, delta.split=FALSE, center.X=TRUE)
REMA.ds <- metahdep(HGU.prep.list, method="REMA", genelist=gene.subset, delta.split=TRUE, center.X=TRUE)


###################################################
### code chunk number 12: metahdep.Rnw:375-377
###################################################
HBLM <- metahdep(HGU.prep.list, method="HBLM", genelist=gene.subset, delta.split=FALSE, center.X=TRUE, two.sided=TRUE)
HBLM.ds <- metahdep(HGU.prep.list, method="HBLM", genelist=gene.subset, delta.split=TRUE, center.X=TRUE, two.sided=TRUE)


###################################################
### code chunk number 13: metahdep.Rnw:395-398
###################################################
pvals <- data.frame(FEMA=FEMA$Intercept.pval, REMA=REMA$Intercept.pval, REMA.ds=REMA.ds$Intercept.pval,
  HBLM=HBLM$Intercept.pval, HBLM.ds=HBLM.ds$Intercept.pval)
plot(pvals, main='Comparison of Intercept P-values from Meta-Analysis',pch=16,cex=.5)


###################################################
### code chunk number 14: metahdep.Rnw:425-428
###################################################
library(metahdep)
data(gloss)
gloss.Table1


###################################################
### code chunk number 15: metahdep.Rnw:435-436
###################################################
gloss.Table1


###################################################
### code chunk number 16: metahdep.Rnw:448-452
###################################################
attach(gloss.Table1)
sp <- sqrt( ((n1-1)*s1^2+(n2-1)*s2^2) / (n1+n2-2) )
c.correct <- (1-3/(4*dfE-1))
theta <- c.correct * (y2-y1)/sp


###################################################
### code chunk number 17: metahdep.Rnw:460-463
###################################################
pE <- c.correct^2*dfE/(dfE-2)
var.theta <- pE*(1/n2+1/n1) + (pE-1)*theta^2
V <- diag(var.theta)


###################################################
### code chunk number 18: metahdep.Rnw:475-480
###################################################
V[2,3] <- V[3,2] <- (pE[2]-1)*theta[2]*theta[3]
V[4,5] <- V[5,4] <- (pE[4]-1)*theta[4]*theta[5]
V[10,11] <- V[11,10] <- (pE[10]-1)*theta[10]*theta[11]
V[10,12] <- V[12,10] <- (pE[10]-1)*theta[10]*theta[12]
V[11,12] <- V[12,11] <- (pE[11]-1)*theta[11]*theta[12]


###################################################
### code chunk number 19: metahdep.Rnw:488-496
###################################################
X <- matrix(0,nrow=18,ncol=6)
X[,1] <- 1
X[Rec==2,2] <- 1
X[Time==2,3] <- 1
X[Part==2,4] <- 1
X[Part==3,5] <- 1
X[Pct==2,6] <- 1
colnames(X) <- c('Intercept','Rec2','Time2','Part2','Part3','Pct2')


###################################################
### code chunk number 20: metahdep.Rnw:518-534
###################################################
# define dependence groups for delta-splitting models
dep.groups <- list( c(2,3,4,5), c(10,11,12) )

# fixed effects meta-analysis
FEMA.indep <- metahdep.FEMA(theta, diag(diag(V)), X, center.X=TRUE)
FEMA.dep <- metahdep.FEMA(theta, V, X, center.X=TRUE)

# random effects meta-analysis
REMA.indep <- metahdep.REMA(theta, diag(diag(V)), X, center.X=TRUE)
REMA.dep <- metahdep.REMA(theta, V, X, center.X=TRUE)
REMA.ds <- metahdep.REMA(theta, V, X, center.X=TRUE, delta.split=TRUE, dep.groups=dep.groups)

# hierarchical Bayes meta-analysis
HBLM.indep <- metahdep.HBLM(theta, diag(diag(V)), X, center.X=TRUE, two.sided=TRUE)
HBLM.dep <- metahdep.HBLM(theta, V, X, center.X=TRUE, two.sided=TRUE)
HBLM.ds <- metahdep.HBLM(theta, V, X, center.X=TRUE, two.sided=TRUE, delta.split=TRUE, dep.groups=dep.groups, n=20, m=20)


###################################################
### code chunk number 21: metahdep.Rnw:548-563 (eval = FALSE)
###################################################
## r1 <- round(c(FEMA.indep$beta.hats[1], FEMA.indep$beta.hat.p.values[1],NA,NA),4)
## r2 <- round(c(FEMA.dep$beta.hats[1], FEMA.dep$beta.hat.p.values[1],NA,NA),4)
## r3 <- round(c(REMA.indep$beta.hats[1], REMA.indep$beta.hat.p.values[1], REMA.indep$tau2.hat, NA),4)
## r4 <- round(c(REMA.dep$beta.hats[1], REMA.dep$beta.hat.p.values[1], REMA.dep$tau2.hat, NA),4)
## r5 <- round(c(REMA.ds$beta.hats[1], REMA.ds$beta.hat.p.values[1], REMA.ds$tau2.hat, REMA.ds$varsigma.hat),4)
## #r5 <- c(0.5261, 0.0026, 0.2936, -0.0816)
## r6 <- round(c(HBLM.indep$beta.hats[1], HBLM.indep$beta.hat.p.values[1],  HBLM.indep$tau.var + HBLM.indep$tau.hat^2, NA),4)
## r7 <- round(c(HBLM.dep$beta.hats[1], HBLM.dep$beta.hat.p.values[1],  HBLM.dep$tau.var + HBLM.dep$tau.hat^2, NA),4)
## r8 <- round(c(HBLM.ds$beta.hats[1], HBLM.ds$beta.hat.p.values[1],  HBLM.ds$tau.var + HBLM.ds$tau.hat^2, HBLM.ds$varsigma.hat),4)
## model <- c(rep('Fixed',2),rep('Random',3),rep('Bayes',3))
## dep <- c('Independent','Dependent',rep(c('Independent','Dependent','Delta-split'),2))
## Table2 <- as.data.frame(cbind(model,dep,rbind(r1,r2,r3,r4,r5,r6,r7,r8)))
## colnames(Table2) <- c('Model','Structure','Intercept','P','tau-square','varsigma')
## rownames(Table2) <- NULL
## Table2


###################################################
### code chunk number 22: metahdep.Rnw:570-585
###################################################
r1 <- round(c(FEMA.indep$beta.hats[1], FEMA.indep$beta.hat.p.values[1],NA,NA),4)
r2 <- round(c(FEMA.dep$beta.hats[1], FEMA.dep$beta.hat.p.values[1],NA,NA),4)
r3 <- round(c(REMA.indep$beta.hats[1], REMA.indep$beta.hat.p.values[1], REMA.indep$tau2.hat, NA),4)
r4 <- round(c(REMA.dep$beta.hats[1], REMA.dep$beta.hat.p.values[1], REMA.dep$tau2.hat, NA),4)
#r5 <- round(c(REMA.ds$beta.hats[1], REMA.ds$beta.hat.p.values[1], REMA.ds$tau2.hat, REMA.ds$varsigma.hat),4)
r5 <- c(0.5261, 0.0026, 0.2936, -0.0816)
r6 <- round(c(HBLM.indep$beta.hats[1], HBLM.indep$beta.hat.p.values[1],  HBLM.indep$tau.var + HBLM.indep$tau.hat^2, NA),4)
r7 <- round(c(HBLM.dep$beta.hats[1], HBLM.dep$beta.hat.p.values[1],  HBLM.dep$tau.var + HBLM.dep$tau.hat^2, NA),4)
r8 <- round(c(HBLM.ds$beta.hats[1], HBLM.ds$beta.hat.p.values[1],  HBLM.ds$tau.var + HBLM.ds$tau.hat^2, HBLM.ds$varsigma.hat),4)
model <- c(rep('Fixed',2),rep('Random',3),rep('Bayes',3))
dep <- c('Independent','Dependent',rep(c('Independent','Dependent','Delta-split'),2))
Table2 <- as.data.frame(cbind(model,dep,rbind(r1,r2,r3,r4,r5,r6,r7,r8)))
colnames(Table2) <- c('Model','Structure','Intercept','P','tau-square','varsigma')
rownames(Table2) <- NULL
Table2