inst/doc/MultiMed.R

In MultiMed: Testing multiple biological mediators simultaneously

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

###################################################
### code chunk number 1: setup
###################################################
library(MultiMed)


###################################################
### code chunk number 2: simSingMed
###################################################
set.seed(20183)

alpha <- 0.2
beta <- 0.2
gamma <- 0.4
n <- 100
sigma2E <- 1
sigma2M <- 1 - alpha^2
sigma2Y <- 1 - beta^2 * (1 - alpha^2) - (alpha * beta + gamma)^2

## exposure:
E <- rnorm(n, 0, sd = sqrt(sigma2E))
## mediator:
M <- matrix(0, nrow = n, ncol = 1)
M[, 1] <- rnorm(n, alpha * E, sd = sqrt(sigma2M))
## outcome:
Y <- rep(0, n)
for (subj in 1:n) Y[subj] <- rnorm(1, beta * M[subj, ], sd = sqrt(sigma2Y))


###################################################
### code chunk number 3: testSingMed
###################################################
medTest(E, M, Y, nperm = 500)


###################################################
### code chunk number 4: simMultMed
###################################################
set.seed(380184)

alpha <- c(rep(0, 6), rep(0.3, 2), rep(0, 2))
beta <- c(rep(0, 6), rep(0, 2), rep(0.3, 2))
gamma <- 0.6

alpha
beta

n <- 100

sigma2E <- 1
sigma2M <- 1-alpha^2
sigma2Y <- 1-sum(beta^2*sigma2M)-(sum(alpha*beta)+gamma)^2

sigma2M
sigma2Y


###################################################
### code chunk number 5: simMultMed2
###################################################
K <- length(alpha)

## exposure:
E <- rnorm(n, 0, sd = sqrt(sigma2E))
## mediator:
M <- matrix(0, nrow = n, ncol = K)
for (i in 1:K) {
  M[, i] <- rnorm(n, alpha[i] * E, sd = sqrt(sigma2M[i]))
}
## outcome:
Y <- rep(0, n)
for (subj in 1:n)
  Y[subj] <- rnorm(1, sum(beta*M[subj,])+gamma*E[subj], sd=sqrt(sigma2Y))


###################################################
### code chunk number 6: testMultMed
###################################################
medTest(E, M, Y, nperm = 500)


###################################################
### code chunk number 7: loadData
###################################################
data(NavyAdenoma)


###################################################
### code chunk number 8: exploreData1
###################################################
colnames(NavyAdenoma)[1:5]


###################################################
### code chunk number 9: exploreData2
###################################################
colnames(NavyAdenoma)[c(6:9,154)]
colnames(NavyAdenoma)[155]
table(NavyAdenoma$Adenoma)


###################################################
### code chunk number 10: getWeights
###################################################
prev <- 0.228

p <- sum(NavyAdenoma$Adenoma==1)/nrow(NavyAdenoma)
p

w <- rep(NA, nrow(NavyAdenoma))
w[NavyAdenoma$Adenoma == 1] <- prev/p
w[NavyAdenoma$Adenoma == 0] <- (1-prev)/(1-p)

table(w)


###################################################
### code chunk number 11: medFish
###################################################
set.seed(840218)
medsFish <- medTest(E=NavyAdenoma$Fish,
                    M=NavyAdenoma[, 6:154],
                    Y=NavyAdenoma$Adenoma,
                    Z=NavyAdenoma[, 2:5],
                    nperm=1000, w=w,
                    useWeightsZ=FALSE)


###################################################
### code chunk number 12: findTopMeds
###################################################
rownames(medsFish) <- colnames(NavyAdenoma[,-c(1:5, 154)])
medsFish[which.min(medsFish[,"p"]),,drop=FALSE]