doc/Summary.R

In wal615/prime.total.effect: Total effect estimation of the environmental mixtures

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "",
  cache = TRUE,
  eval = TRUE
)

## ----GCTA.rr, echo = T, cache=TRUE--------------------------------------------
library(prime.total.effect)
# Prepare data
data <- hemoglobin.PCB[,1:39] # only select y and PCBs 
data <- na.omit(data) # remove all the NA 
x <- std.fn(data[,-1])   # standardize PCBs
y <- data$LBXGH

# Variance of y

var(y)

# Calculate the total effect by GCTA method

res.GCTA.beta2 <- GCTA.rr(y = y, x = x, target = "beta2", bs.iter = 1, verbose = FALSE) # "beta2" means total effect
res.GCTA.beta2$beta2
res.GCTA.r2 <- GCTA.rr(y = y, x = x, target = "r2", bs.iter = 1, verbose = FALSE) # r2 means precentage
res.GCTA.r2$r2

# Add the bootstrap 
# bs.iter means the number of bootstrap sample. 
# In example, we only use 10 iterations but for real case it should be larger than 1000
res.GCTA.beta2 <- GCTA.rr(y = y, x = x, target = "beta2", bs.iter = 10, verbose = TRUE) # "beta2" means total effect
res.GCTA.beta2
res.GCTA.r2 <- GCTA.rr(y = y, x = x, target = "r2", bs.iter = 10, verbose = TRUE) # r2 means precentage
res.GCTA.r2

## ----decorrelation, cache=TRUE------------------------------------------------
library(prime.total.effect)
# Prepare data
data <- hemoglobin.PCB[,1:39]         # only select y and PCBs 
data <- na.omit(data)                 # remove all the NA 
x <- std.fn(data[,-1])   # standardize PCBs
y <- data$LBXGH

# decorrelation step
x.hist <- empirical.cov.method(x)$uncorr.data # decorrelation by directly by the sample covariance matrix when n >p 
x.glasso <- GLASSO.method(x)$uncorr.data # decorrelation by glasso method. note it could also be used when n < p 

cor(x)[1:5,1:5]
cor(x.hist)[1:5,1:5]
cor(x.glasso)[1:5,1:5]

# compare the heritability estimations among different decorrelation method

res.GCTA.r2 <- GCTA.rr(y = y, x = x, target = "r2", bs.iter = 1, verbose = FALSE) # r2 means precentage
res.GCTA.r2$r2
res.GCTA.r2.hist <- GCTA.rr(y = y, x = x.hist, target = "r2", bs.iter = 1, verbose = FALSE) # r2 means precentage
res.GCTA.r2.hist$r2
res.GCTA.r2.glasso <- GCTA.rr(y = y, x = x.glasso, target = "r2", bs.iter = 1, verbose = FALSE) # r2 means precentage
res.GCTA.r2.glasso$r2

## ----add interaction terms, cache=TRUE----------------------------------------
library(prime.total.effect)
# Prepare data
data <- hemoglobin.PCB[,1:39]         # only select y and PCBs 
data <- na.omit(data)                 # remove all the NA 
x <- std.fn(data[,-1])   # standardize PCBs
y <- data$LBXGH

# add two-way interaction terms
x.with.inter <- std.fn(add.inter(x))
colnames(x.with.inter)[c(1:5,38:43)]

# Decorrelation step
x.hist <- empirical.cov.method(x.with.inter)$uncorr.data # decorrelation by directly by the sample covariance matrix

# Calculate the precentage
res.GCTA.r2.hist <- GCTA.rr(y = y, x = x.hist, target = "r2", bs.iter = 1, verbose = FALSE) # r2 means precentage
res.GCTA.r2.hist$r2

## ----master function example 1, cache=TRUE------------------------------------
library(prime.total.effect)

# Prepare data
data <- hemoglobin.PCB[,1:39]         # only select y and PCBs 
data <- na.omit(data)                 # remove all the NA 
x <- std.fn(data[,-1])   # standardize PCBs
y <- data$LBXGH

# select the GCTA methods
total.effect.method <- GCTA.rr
# total.effect.args.list contains all the arguments for this method, the details of arguments can be found in `help(GCTA.rr)`
total.effect.args.list <- list(target = "r2", bs.iter = 10, verbose = FALSE)

total.effect.test.res <- total.effect.master(x = x, y = y, inter = TRUE,
                                             decorrelation.method = "historical",
                                             total.effect.method = total.effect.method, 
                                             total.effect.args.list = total.effect.args.list)
total.effect.test.res

## ----master function example 2, cache=TRUE------------------------------------
# subset 200 observation for the original data set
# so that after adding the interaction terms, we have the n < p case 
# then we could use the EigenPrism and Dickers method 
data <- data[1:500, ]
x <- std.fn(data[,-1])   # standardize PCBs
y <- data$LBXGH

# select the EigenPrism methods
total.effect.method <- EigenPrism 
total.effect.args.list <- list(target = "r2", alpha = 0.05) # the details of the arguments are in `EigenPrism`

total.effect.test.res <- total.effect.master(x = x, y = y, inter = TRUE,
                                             decorrelation.method = "glasso", rho = 0.01, thr = 0.01,
                                             total.effect.method = total.effect.method, 
                                             total.effect.args.list = total.effect.args.list)
total.effect.test.res

# select the Dicker methods
total.effect.method <- Dicker.2013 
total.effect.args.list <- list(target = "r2") # the details of the arguments are in `EigenPrism`

total.effect.test.res <- total.effect.master(x = x, y = y, inter = TRUE,
                                             decorrelation.method = "glasso", rho = 0.01, thr = 0.01,
                                             total.effect.method = total.effect.method, 
                                             total.effect.args.list = total.effect.args.list)
total.effect.test.res