Nothing
R/wga_GbyE.R
In CGEN: An R package for analysis of case-control studies in genetic epidemiology
Defines functions
meta.GxE
getCov12Object
getCovObject
getBetaObject
meta.EB
get.EB.A_matrix
meta.fixed
scan.UML_CML
extractEst
extract_UML.CML
outNames.cov2
outNames.cov
outNames.me
writeOut
UML_CML_GxE_parms
apply_zero.vars
logistic.dsgnMat
logistic.vnames
snp.main
snp.logistic
snp.scan.logistic
sMatrix.logistic
Documented in
snp.logistic
# History Apr 21 2008 Initial coding
# May 09 2008 Use returned hessian to compute cov
# May 29 2008 Redo the log-likelihood calculation
# May 31 2008 Move tests to GbyE function.
# Jun 02 2008 Add gradient function for the optimizer
# Add code for new output file.
# Jun 03 2008 Add test2.vars for a user specified Wald test.
# Jun 04 2008 Add tests.init (standard logistic reg.)
# Jun 06 2008 Add parameter estimates and standard errors
# to out.file.
# Make vnames a seperate function
# Jun 07 2008 Combine GbyE.1snp and GbyE.irls
# Jun 09 2008 Change input argument names, function names.
# Remove strataMat from input.
# Jun 13 2008 Catch errors when calling logistic.SNP
# Jun 16 2008 Rename getLoglike to getLoglike.glm
# Jun 17 2008 Rename strata vars to Allele_freq
# Put strata parms in a separate object
# Jun 18 2008 Update for change in getWaldTest
# Jun 20 2008 Change X.vars to X.main, V.vars to X.int
# Jun 30 2008 Remove call to setUpSummary
# Jul 02 2008 Add code for streamed input
# Jul 09 2008 Change in getEB function:
# cmat <- cbind(diag(temp), diag(1-temp))
# Jul 11 2008 Rename the functions to snp.logistic and
# snp.scan.logistic
# Jul 11 2008 Add output for snp and subject counts
# Jul 17 2008 Add code for factor snps
# Jul 18 2008 Add test.inter
# Jul 18 2008 Add code for saving individual tests and p-values
# Jul 21 2008 Add option to output corr parms
# Jul 22 2008 Change code to call getTest function
# Jul 26 2008 Add tests and tests.1df options
# The tests option will replace the test2.vars option
# Aug 07 2008 Check for constant variables in getDsgnMat
# Aug 11 2008 Define classes for return objects
# Aug 21 2008 Add option for joint and stratified effects
# Aug 25 2008 Rename snp.logistic to snp.main
# Add snp.logistic function
# Aug 29 2008 Update for change in getWaldTest
# Oct 08 2008 Allow for different family if only.UML = 1
# Oct 10 2008 Allow for no recoding of genotypes when stream = 1
# Oct 21 2008 Make more efficient for only.UML = 1
# Oct 29 2008 Change the line in snp.scan.logistic from
# snp <- as.integer(snp) to
# snp <- as.numeric(snp)
# Dec 10 2008 Allow the user to specify formulas in
# snp.logistic.
# Jan 30 2009 Add option to output UML base model
# Mar 05 2009 Get MAF from getData.1 function (stream = 0)
# Mar 05 2009 Update for stream = 1
# Mar 16 2009 Do not set cc.var to response.var by default
# Jun 16 2009 Print out response levels only for binomial family
# Jul 10 2009 Check if design.V0 is NULL in getInter.Vars()
# Jul 22 2009 Move temp.list inside the options list
# Aug 03 2009 Change in outputRow to match the change in
# effects.init return list
# Sep 17 2009 Add in output option in getDelimiter in snp.scan.logistic
# Change as.integer to as.numeric
# Add errorCheck option in snp.main
# Oct 01 2009 Check variable names in snp.logistic
# Oct 16 2009 Change in recode.geno
# Oct 16 2009 Add major/minor allele to out file in snp.scan.logistic
# Oct 20 2009 Fix bug in getting response0 in snp.scan.logistic
# Sometimes it is character
# Oct 28 2009 Make sure design matrices from logistic.dsgnMat
# are numeric
# Dec 17 2009 Remove IRLS algorithm
# Add code for genetic.model = 3
# Dec 23 2009 Generalize the stratification
# Jan 04 2010 Update snp.logistic for missing values
# Feb 01 2010 Wrap callOptim in try function in snp.main
# Feb 04 2010 In snp.logistic, check for 1 strata.var that is constant
# Feb 12 2010 Use chol and chol2inv instead of solve.
# Add option for C code.
# Feb 17 2010 Check CML variances from returned C code
# Feb 17 2010 Make getInit function in snp.main more efficient
# Feb 26 2010 Use the snp variable name in snp.logistic
# Mar 04 2010 Use normVarNames function in snp.logistic
# Mar 14 2010 Add option for using (2) only controls to determine
# the major/minor allele or (1) all subjects
# The default is 1
# Mar 18 2010 Add package="CGEN" option
# Mar 29 2010 Add return values to snp.logistic
# Mar 30 2010 Add option for fixing parameters
# Apr 02 2010 Add option zero.vars to snp.logistic
# Apr 05 2010 Change call to C code
# Apr 14 2010 Change default reltol to 1e-8
# Apr 22 2010 Fix bug in getInit with adding snp back when zeroSNP = 1
# Apr 26 2010 Change code for when there are NAs with UML parms
# Nov 04 2011 Check if snp.var is a main effect
# Dec 23 2011 Add function in snp.main to check CML initial estimates
# Allow for initial estimates to be passed in
# Jul 18 2012 Remove errorCheck option in snp.main. Was causing a problem
# for different genetic models.
# Feb 02 2013 Add warning message for continuous variables
# Mar 26 2013 Add code for imputed SNPs
# Mar 29 2013 Add new functions for outputing UML-CML estimates
# Apr 03 2013 Add code for meta-analysis
# Apr 04 2013 Generalize code for UML-CML estimates
# Apr 05 2013 Modify snp.scan.logistic for imputed snps
# Sep 01 2013 Allow snp to be a 3 column matrix in snp.main
# Oct 22 2013 Return full CML covariance matrix from snp.main
# Oct 25 2013 Add function to return strat matrix
######################################################################
# TO DO:
# Modify the different tests for the different genetic models
# Update test2 code for factors
# Check for consistent input arguments (out.est, only.UML,...)
######################################################################
sMatrix.logistic <- function(data, strata.var, facVars) {
sflag <- !is.null(strata.var)
sform <- "formula" %in% class(strata.var)
nobs <- nrow(data)
svarCat1 <- FALSE
# Check for constant strata variable
if ((sflag) && (!sform) && (length(strata.var) == 1)) {
svarCat1.v <- data[, strata.var]
svarCat1.s <- sort(unique(svarCat1.v))
svarCat1.n <- length(svarCat1.s)
if (svarCat1.n == 1) sflag <- FALSE
if (sflag) {
if ((is.character(svarCat1.v)) || (is.factor(svarCat1.v))) svarCat1 <- TRUE
}
}
if (!sflag) {
design.S0 <- matrix(data=1, nrow=nobs, ncol=1)
s.1catVar <- 1
} else {
if (svarCat1) {
design.S0 <- matrix(data=0, nrow=nobs, ncol=svarCat1.n)
s.1catVar <- 1
for (i in 1:svarCat1.n) design.S0[(svarCat1.v %in% svarCat1.s[i]), i] <- 1
} else {
design.S0 <- logistic.dsgnMat(data, strata.var, facVars, removeInt=0)$designMatrix
s.1catVar <- 0
}
}
list(design.S0=design.S0, s.1catVar=s.1catVar)
} # END: sMatrix.logistic
# Function to perform SNP by environment interaction analysis.
snp.scan.logistic <- function(snp.list, pheno.list, op=NULL) {
# INPUT:
# snp.list (See snp.list.wordpad)
#################################################################
# pheno.list (See pheno.list.wordpad)
# This list must include the name "response.var",
# and can also include the names "strata.var",
# "main.vars", "int.vars" (see below)
# response.var Name of the response variable. This variable
# must be coded as 0 and 1.
# No default
# strata.var Stratification variable name(s) or a formula.
# The default is NULL so that all observations
# belong to the same strata.
# main.vars Vector of variables names in the file pheno.list$file
# that will be included in the model as main effects.
# main.vars can also be a formula.
# The default is NULL.
# int.vars Vector of variable names in the file pheno.list$file
# that will interact with each snp in the snp data.
# int.vars can also be a formula.
# The default is NULL.
# cc.var The name of the case-control variable. This variable must
# be coded as 0-1. The variable is used to determine
# the MAF based on controls.
# The default is NULL.
#################################################################
# op List with the following names.
# genetic.model 0-2
# 0: trend
# 1: dominant
# 2: recessive
# 3: general
# The default is 0.
# reltol Stopping tolerance
# The default is 1e-8
# maxiter Maximum number of iterations
# The default is 100
# optimizer One of : "BFGS", "Nelder-Mead", "BFGS", "CG",
# "L-BFGS-B", "SANN".
# The default is "BFGS"
# print 0 or 1 to print each list of output
# The default is 0
# test.omnibus 0 or 1 for the test
# involving the snp and interaction terms.
# The default is 1
# test.main 0 or 1 for computing
# the test for the main effect of the snp.
# The default is 1.
# test.inter 0 or 1 for computing
# the test for the interaction terms.
# The default is 0.
# tests List of character vector of variable names that
# will be used in Wald tests.
# The default is NULL.
# tests.1df Character vector of variable names to use for 1 df
# Wald tests.
# The default is NULL.
# tests.UML NULL or a list with at least one of the following
# names: "omnibus", "main", or "inter". These names
# can be set to one of the following values:
# "WALD", "LR", or "SCORE".
# Example: test.UML=list(main="WALD", omnibus="LR") will
# compute a wald test for the main effect of SNP and
# a likelihood ratio test for the SNP and interactions.
# !!!! Only for family=binomial !!!!!!
# The default value for each name is "WALD"
# The default is NULL.
# tests.CML 0 or 1 for computing the WALD tests using the
# conditional ML estimates.
# The default is 1.
# tests.EB 0 or 1 for computing the WALD tests using the
# empirical bayes estimates
# The default is 0.
# geno.counts 0 or 1 to add the genotype counts to out.file
# The default is 0.
# subject.counts 0 or 1 to add the number of cases and controls to
# out.file.
# The default is 0.
# allele.cc 1 or 2 to use all subjects
###############################################################
# effects List for joint/stratified effects
# Names in the list must be "var", "type",
# "var.levels", "snp.levels", "var.base"
# The default is NULL.
# var Variable name to compute the effects with the SNP
# No default
# type 1 or 2 or c(1, 2) , 1 = joint, 2 = stratified
# The default is 1.
# var.levels (Only for continuous var) Numeric vector of the
# levels to be used in the calculation.
# The default is 0.
# var.base (Only for continuous var) Baseline level.
# The default is 0
# snp.levels One of the following: 0, 1, 2, c(0,1), c(0,2),
# c(1,2) or c(0, 1, 2)
# The default is 1.
# method Character vector containing any of the following:
# "UML", "CML", "EB"
# The default is c("UML", "CML", "EB")
###############################################################
# out.est List with the names "parms" and "method".
# Use this option if you want to save parameter
# estimates and standard errors to out.file.
# parms 1-3 or a character vector of the parameters to
# save statistics on.
# 1: Main effect is saved.
# 2: Main effect and interactions are saved.
# 3: All parameters are saved.
# method Character vector containing any of the following:
# "UML", "CML", "EB"
# For example, setting parms=3 and method=c("CML", "EB")
# will write all the parameter estimates and
# standard errors for the methods CML and EB.
# The default value of out.est is NULL.
# what Character vector containing any of the following:
# "beta", "se", "test", "pvalue"
# corr.parms List of character vectors of length 2 giving the
# pairs of variables for correlations.
###############################################################
# out.file NULL or file name to save summary information for
# each snp. The output will at least contain the columns
# "SNP" and "MAF". MAF is the minor allele frequency
# form the controls. Additional columns in this file
# are based on the values of test.omnibus, test.main,
# test2.vars, tests.UML, tests.CML, tests.EB, and
# out.est.
# The default is NULL.
# out.dir NULL or the output directory to store the output
# lists for each SNP. A seperate file will be created
# for each snp in the snp data set.
# The file names will be the out_<snp>.rda
# The object names are called "ret".
# The default is NULL.
# base.outfile NULL or path to base model rda file
# The default is NULL.
#####################################################################
# temp.list See temp.list.doc
#####################################################################
# RETURN:
# The list returned is from the last analysis performed.
# Function to call different test functions
getTest <- function(fit, vars, test, which, model) {
# fit List of parms, cov, and loglike
# vars Variables to test
# test "WALD", "LR", or "SCORE"
# which 1-4 1=omnibus, 2=main, 3=inter, 4=test2
# model "UML", "CML", or "EB"
if (test == "WALD") {
# Set up a list for calling getWaldTest
temp <- list(parms=fit$parms, cov=fit$cov)
return(getWaldTest(temp, vars))
}
# Determine if test is valid
if (test == "LR") {
# LR is not for EB
if (model == "EB") return(list(test=NA, df=NA, pvalue=NA))
if ((model == "CML") && (which %in% c(1, 2))) return(list(test=NA, df=NA, pvalue=NA))
} else if (test == "SCORE") {
if ((model == "CML") && (which %in% c(1, 2))) return(list(test=NA, df=NA, pvalue=NA))
}
# Determine if interactions are in the model
if (which == 2) {
X.int <- design.V
xflag <- 1
} else {
X.int <- NULL
xflag <- 0
}
# Factor snp if genetic.model = 3
if (genetic.model == 3) snp <- factor(snp)
# Determine if SNP is in the model
if (which == 3) {
int.vec <- snp
flag <- 1
} else {
int.vec <- NULL
flag <- 0
}
# Base model
if (model == "UML") {
fit0 <- callGLM(response, X.main=design.X, X.int=X.int, int.vec=snp,
family="binomial", prefix="SNP_", retX=TRUE,
retY=TRUE, inc.int.vec=flag)
} else {
fit0 <- snp.main(response, snp, X.main=design.X, X.int=design.V,
X.strata=design.S, op=op)
fit0 <- fit0[[model]]
# Add info for score tests
if (test == "SCORE") {
# Add x, y, linear.predictors
temp <- getXBeta(cbind(design.X, int.vec, X.int), fit0$parms)
fit0$linear.predictors <- temp$XBeta
fit0$x <- temp$X
fit0$y <- response
} else {
# LR test, add the rank
fit0$rank <- ncol(fit0$cov)
}
}
if (test == "LR") {
# Likelihood ratio test. Add the rank
fit$rank <- ncol(fit$cov)
return(likelihoodRatio(fit, fit0))
} else {
# Score test
# Get the matrix to test
temp <- NULL
if ((genetic.model == 3) && (!xflag || !flag)) {
snp <- createDummy(snp)$data
}
if (!xflag) temp <- addInterVars(NULL, snp, design.V)$data
if (!flag) temp <- cbind(temp, snp)
return(score.logReg(fit0, temp))
}
} # END: getTest
# Function to add tests to the return list
addToList <- function(ret, which) {
# ret Return list
# which "UML", "CML", "EB"
if (which == "UML") {
# tests.UML.vec is actually a list
test <- tests.UML.vec
} else {
test <- rep("WALD", times=3)
}
field <- paste(which, c(".omnibus", ".main", ".inter"), sep="")
temp <- ret[[which]]
if (!is.null(temp)) {
if (omniFlag) ret[[field[1]]] <-
getTest(temp, omni.vars[[vListIndex]], test[1], 1, which)
if (mainFlag) ret[[field[2]]] <-
getTest(temp, main.vars[[vListIndex]], test[2], 2, which)
if (interFlag) ret[[field[3]]] <-
getTest(temp, inter.vars[[vListIndex]], test[3], 3, which)
# tests will start from field 4
if (test2Flag) {
field2 <- paste(which, ".test", 1:n.tests, sep="")
for (i in 1:n.tests) {
ret[[field2[i]]] <- getTest(temp, test2.vars[[i]], "WALD", 4, which)
}
}
} else {
if (omniFlag) ret[[field[1]]] <- list(test=NA, pvalue=NA, df=NA)
if (mainFlag) ret[[field[2]]] <- list(test=NA, pvalue=NA, df=NA)
if (interFlag) ret[[field[3]]] <- list(test=NA, pvalue=NA, df=NA)
if (test2Flag) {
temp <- list(test=NA, pvalue=NA, df=NA)
for (i in 1:n.tests) ret[[field2[i]]] <- temp
}
}
ret
} # END: addTest
# Function to tests to the output vector
addToVec <- function(outVec, which) {
# outVec Output vector
# which "omnibus", "main", "inter", or "testi"
tname <- paste(which, ".test", sep="")
pname <- paste(which, ".pvalue", sep="")
dname <- paste(which, ".df", sep="")
if (tests.UML) {
t <- paste("UML.", tname, sep="")
p <- paste("UML.", pname, sep="")
d <- paste("UML.", dname, sep="")
r <- paste("UML.", which, sep="")
outVec[t] <- ret[[r]]$test
outVec[p] <- ret[[r]]$pvalue
outVec[d] <- ret[[r]]$df
}
if (tests.CML) {
t <- paste("CML.", tname, sep="")
p <- paste("CML.", pname, sep="")
d <- paste("CML.", dname, sep="")
r <- paste("CML.", which, sep="")
outVec[t] <- ret[[r]]$test
outVec[p] <- ret[[r]]$pvalue
outVec[d] <- ret[[r]]$df
}
if (tests.EB) {
t <- paste("EB.", tname, sep="")
p <- paste("EB.", pname, sep="")
d <- paste("EB.", dname, sep="")
r <- paste("EB.", which, sep="")
outVec[t] <- ret[[r]]$test
outVec[p] <- ret[[r]]$pvalue
outVec[d] <- ret[[r]]$df
}
outVec
} # END: addToVec
# Function for the number of cases and controls
getCCcounts <- function(notMiss) {
# Note: The cntrl vector equals 1 for controls.
# So the returned vector from table(), will have the counts
# for the cases and then the controls, which is the order in
# the output file.
if (is.null(notMiss)) return(subj.cnts0)
temp <- cntrl[notMiss]
nn <- sum(temp)
return(c(length(temp)-nn, nn))
} # END: getCCcounts
# Function to get variables for main test
getMain.vars <- function() {
v2 <- NULL
v1 <- getVarNames.snp(prefix=op$snpName, genetic.model=genetic.model)
if (genetic.model == 3) {
v2 <- getVarNames.snp(prefix=op$snpName, genetic.model=0)
}
list(v1, v2)
} # END: getMain.vars
# Function to get variables for omnibus test
getOmni.vars <- function() {
# Combine main effect and interaction vars
mvars <- getMain.vars()
temp <- getInter.vars()
v1 <- c(mvars[[1]], temp[[1]])
v2 <- c(mvars[[2]], temp[[2]])
list(v1, v2)
} # END: getMain.vars
# Function to get the interaction vars
getInter.vars <- function() {
if (is.null(design.V0)) return(NULL)
v2 <- NULL
v1 <- getVarNames.int(design.V0, prefix=op$snpName, genetic.model=genetic.model, sep=":")
if (genetic.model == 3) {
v2 <- getVarNames.int(design.V0, prefix=op$snpName, genetic.model=0, sep=":")
}
list(v1, v2)
} # END: getInter.vars
# Function to return the test2 vars
getTest2.vars <- function() {
tests <- getListName(op, "tests")
ret <- list()
# Loop over each set of variables
for (i in 1:length(tests)) {
testi <- tests[[i]]
# Determine if any of the test2 variables are factors
temp <- 0
if (facFlag) temp <- testi %in% facVars
if (!any(temp)) {
test2.vars <- testi
} else {
# First add the variables that are not factors
test2.vars <- testi[as.logical(1-temp)]
# Now add the dummy variables for the factors
temp <- testi[temp]
for (temp2 in temp) {
# Get the dummy variable names in one of the lists
temp3 <- getListName(X.newVars, temp2)
if (!is.null(temp3)) {
test2.vars <- c(test2.vars, temp3)
} else {
temp3 <- getListName(V.newVars, temp2)
test2.vars <- c(test2.vars, temp3)
}
}
}
# Check the variable names
temp <- check.vec.char(test2.vars, "tests", len=NULL,
checkList=log.vnames$all)
if (temp) stop("ERROR with op$tests")
# Add to the list
ret[[i]] <- test2.vars
}
ret
} # END: getTest2.vars
# Function to return the index to use to get the test variable names
getVListIndex <- function(n) {
if (genetic.model != 3) {
return(1)
} else {
if (n == 3) {
return(1)
} else {
return(2)
}
}
} # END: getVListIndex
# Function to in initialize the output file
initOutfile <- function() {
# Create an output vector
vnames <- NULL
vec <- c(".test", ".pvalue", ".df")
if (omniFlag) {
temp <- paste("omnibus", vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
if (mainFlag) {
temp <- paste("main", vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
if (interFlag) {
temp <- paste("inter", vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
if (test2Flag) {
for (i in 1:n.tests) {
temp <- paste("test", i, vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
}
est.se <- NULL
est.test <- NULL
est.pval <- NULL
est.corr <- NULL
if (out.est.flag) {
if ((out.est$parms == 2) && (!Vflag)) out.est$parms <- 1
# Vector to hold parm names
what <- out.est$what
temp <- NULL
if (out.est.beta) temp <- est.p[[1]]
if (out.est.se) {
est.se <- list()
for (i in 1:length(est.parms)) {
est.se[[i]] <- paste(est.parms[[i]], ".se", sep="")
}
temp <- c(temp, est.se[[1]])
}
if (out.est.test) {
est.test <- list()
for (i in 1:length(est.parms)) {
est.test[[i]] <- paste(est.parms[[i]], ".test", sep="")
}
temp <- c(temp, est.test[[1]])
}
if (out.est.pval) {
est.pval <- list()
for (i in 1:length(est.parms)) {
est.pval[[i]] <- paste(est.parms[[i]], ".pvalue", sep="")
}
temp <- c(temp, est.pval[[1]])
}
if (out.est.corr) {
est.corr <- NULL
corrs <- getListName(out.est, "corr.parms")
for (i in 1:length(corrs)) {
var <- paste(corrs[[i]], collapse="_", sep="")
est.corr <- c(est.corr, var)
}
temp <- c(temp, est.corr)
}
for (method in out.est$method) {
temp2 <- paste(method, ".", temp, sep="")
vnames <- c(vnames, temp2)
}
}
# Effects
if (effectsFlag) {
names <- effects$out.names
nr <- nrow(names)
for (method in effects$method) {
for (type in effects$type) {
for (i in 1:nr) {
temp <- paste(method, ".eff", type, ".",
names[i,], sep="")
vnames <- c(vnames, temp)
}
# Standard errors
for (i in 1:nr) {
temp <- paste(method, ".eff", type, ".",
names[i,], ".se", sep="")
vnames <- c(vnames, temp)
}
}
}
}
outVec <- double(length(vnames))
names(outVec) <- vnames
# geno and subject counts
if (subj.counts) vnames <- c("N.cases", "N.controls", vnames)
if (geno.counts) vnames <- c("Hom.common", "Heter", "Hom.uncommon", vnames)
vnames <- c("SNP", "Alleles", "MAF", vnames)
fid <- writeVecToFile(vnames, op$out.file, colnames=NULL, type=3,
close=0, sep="\t")
list(fid=fid, outVec=outVec, est.se=est.se,
est.test=est.test, est.pval=est.pval, est.corr=est.corr)
} # END: initOutfile
# Function to create a design matrix
getDsgnMat <- function(vars) {
design <- removeOrKeepCols(phenoData0, vars, which=1)
newVars <- NULL
if (facFlag) {
temp <- vars %in% facVars
if (any(temp)) {
temp <- vars[temp]
temp <- createDummy(design, vars=temp)
design <- temp$data
newVars <- temp$newVars
}
}
design <- as.matrix(design)
# Check for constant variables
design <- checkForConstantVar(design, msg=1)$data
list(designMatrix=design, newVars=newVars)
} # END: getDsgnMat
# Function to output a row to the output file
outputRow <- function() {
cat(snp.name, majMin, maf, "", file=fid, sep="\t")
if (geno.counts) cat(genoCounts, "", file=fid, sep="\t")
if (subj.counts) cat(getCCcounts(subj.notMiss), "", file=fid, sep="\t")
if (omniFlag) outVec <- addToVec(outVec, "omnibus")
if (mainFlag) outVec <- addToVec(outVec, "main")
if (interFlag) outVec <- addToVec(outVec, "inter")
if (test2Flag) {
for (i in 1:n.tests) {
outVec <- addToVec(outVec, paste("test", i, sep=""))
}
}
if (out.est.flag) {
parms <- est.p[[vListIndex]]
for (method in out.est$method) {
temp <- getListName(ret, method)
if (out.est.beta) temp1 <- paste(method, ".", est.parms[[vListIndex]], sep="")
if (out.est.se) temp2 <- paste(method, ".", est.se[[vListIndex]], sep="")
if (out.est.test) temp3 <- paste(method, ".", est.test[[vListIndex]], sep="")
if (out.est.pval) temp4 <- paste(method, ".", est.pval[[vListIndex]], sep="")
if (out.est.corr) temp5 <- paste(method, ".", est.corr, sep="")
if (!is.null(temp)) {
# NOTE: Names that are not found in a named vector
# get assigned a missing value (NA) to them.
se <- sqrt(diag(temp$cov))[parms]
p <- temp$parms[parms]
if (out.est.beta) outVec[temp1] <- p
if (out.est.se) outVec[temp2] <- se
if (out.est.test) outVec[temp3] <- p/se
if (out.est.pval) outVec[temp4] <- 2*pnorm(abs(p/se), lower.tail=FALSE)
if (out.est.corr) {
# The order is the same as est.corr
corrs <- out.est$corr.parms
cov <- temp$cov
for (i in 1:length(corrs)) {
outVec[temp5[i]] <- cov[corrs[[i]][1], corrs[[i]][2]]
}
}
} else {
if (out.est.beta) outVec[temp1] <- NA
if (out.est.se) outVec[temp2] <- NA
if (out.est.test) outVec[temp3] <- NA
if (out.est.pval) outVec[temp4] <- NA
if (out.est.corr) outVec[temp5] <- NA
}
} # END: for (method in out.est$method)
} # END: if (out.est.flag)
# Effects
if (effectsFlag) {
names <- effects$out.names
nr <- nrow(names)
lnames <- effects$list.names
for (method in effects$method) {
i <- 0
for (type in effects$type) {
i <- i + 1
# Get the list of effects and standard errors
temp <- paste(method, lnames[i], sep="")
leff <- ret[[temp]]
if (is.null(leff)) next
eff <- leff$logEffects
for (j in 1:nr) {
temp <- paste(method, ".eff", type, ".", names[j, ], sep="")
outVec[temp] <- eff[j, ]
}
# Standard errors
eff <- leff$logEffects.se
for (j in 1:nr) {
temp <- paste(method, ".eff", type, ".", names[j, ], ".se", sep="")
outVec[temp] <- eff[j, ]
}
}
}
} # END: if (effectsFlag)
cat(outVec, file=fid, sep="\t")
cat("\n", file=fid)
0
} # END: outputRow
# Function to check the out.est list
check.out.est <- function(out.est) {
corr <- getListName(out.est, "corr.parms")
flag <- !is.null(corr)
new <- corr
if (flag) {
if (class(corr) != "list") {
corr <- unique(corr)
n <- length(corr)
if (n > 1) {
index <- 1
# Convert into a list
new <- list()
for (i in 1:(n-1)) {
for (j in (i+1):n) {
list[[index]] <- c(corr[i], corr[j])
index <- index + 1
}
}
} else {
new <- NULL
}
}
}
out.est$corr.parms <- new
if (!flag) {
out.est <- default.list(out.est, c("parms", "method", "what"),
list(1, "CML", c("beta", "se")),
checkList=list(NA, c("UML", "CML", "EB"),
c("beta", "se", "test", "pvalue")) )
}
out.est
} # END: check.out.est
# Function to check the UML test list
check.tests.uml <- function(u) {
if (is.null(u)) return(NULL)
if (!is.list(u)) {
if (u) {
return(list(omnibus="WALD", main="WALD", inter="WALD"))
} else {
return(NULL)
}
}
tests <- c("WALD", "LR", "SCORE")
omni <- getListName(u, "omnibus")
main <- getListName(u, "main")
inter <- getListName(u, "inter")
temp <- c(omni, main, inter)
temp <- match(temp, tests)
if (any(is.na(temp))) stop("ERROR with tests.UML")
# Get the correct order
if (is.null(omni)) omni <- "WALD"
if (is.null(main)) main <- "WALD"
if (is.null(inter)) inter <- "WALD"
u <- list(omnibus=omni, main=main, inter=inter)
u
} # END: check.tests.uml
# Function to check that the options are consistent
checkOptions <- function(op) {
# Check for the tests.1df option
tests.1df <- getListName(op, "tests.1df")
if (!is.null(tests.1df)) {
op$out.est <- list(parms=tests.1df, what=c("test", "pvalue"))
op$tests.1df <- NULL
}
# Check tests list
tests <- getListName(op, "tests")
if (!is.null(tests)) {
# List with vectors is also of type vector
if (!is.list(tests)) tests <- list(tests)
}
op$tests <- tests
# Make sure at least 1 test is being done
temp <- op$test.omnibus + op$test.main + op$test.inter +
as.numeric(!is.null(getListName(op, "out.est"))) +
as.numeric(!is.null(getListName(op, "tests")))
if (!temp) op$test.omnibus <- 1
op$effects <- checkEffects(getListName(op, "effects"))
op
} # END: checkOptions
# Function to check the effects list
checkEffects <- function(eff) {
if (is.null(eff)) return(NULL)
eff <- default.list(eff,
c("var", "type", "var.levels", "snp.levels", "method", "var.base"),
list("ERROR", 1, 0, 1, c("UML", "CML", "EB"), 0),
error=c(1, 0, 0, 0, 0, 0))
temp <- eff$method %in% c("UML", "CML", "EB")
if (any(temp)) {
eff$method <- eff$method[temp]
} else {
eff$method <- c("UML", "CML", "EB")
}
eff
} # END: checkEffects
# Function to compute the effects
calcEffects <- function(eff, ret) {
for (method in eff$method) {
temp <- ret[[method]]
if (!is.null(temp)) {
for (type in eff$type) {
if (type == 1) {
name <- paste(method, ".JointEffects", sep="")
} else {
name <- paste(method, ".StratEffects", sep="")
}
ret[[name]] <- try(effects.init(temp$parms, temp$cov, eff$var1,
eff$var2, eff$snp.levels, eff$var.levels,
base1=0, base2=eff$var.base, int.var=eff$int.var,
effects=type, sep1=eff$sep), silent=TRUE)
if (class(ret[[name]]) == "try-error") ret[[name]] <- NULL
}
}
}
ret
} # END: calcEffects
# Function to set up effects list
setupEffects <- function(eff, facVars, X.vars, V.vars) {
if (is.null(eff)) return(NULL)
# Call before phenoData0 is removed
eff$var1 <- getMain.vars()[[1]]
temp <- nchar(eff$var1)
if (substring(eff$var1, temp, temp) == "_") {
eff$sep <- ""
} else {
eff$sep <- "_"
}
# Check if the variable is a main effect
temp <- eff$var %in% X.vars
if (!any(temp)) {
temp <- paste(eff$var, " is not a main effect. No effects will be computed",
sep="")
warning(temp)
return(NULL)
}
eff$var <- eff$var[temp]
nvar <- length(eff$var)
if ((nvar == 1) && (eff$var %in% facVars)) {
temp <- levels(factor(phenoData0[, eff$var]))
# For now, first level is the baseline
baseline <- paste(eff$var, "_", temp[1], sep="")
temp <- temp[-1]
# Get the dummy variable names
eff$var2 <- paste(eff$var, "_", temp, sep="")
flag <- 1
} else {
eff$var2 <- eff$var
flag <- 0
}
# Check if var is also an interaction
temp <- eff$var %in% V.vars
if (all(temp)) {
if (!flag) {
# For continuous var
eff$int.var <- paste(eff$var1, eff$var, sep="")
} else {
eff$int.var <- paste(eff$var1, eff$var2, sep="")
}
} else {
eff$int.var <- NULL
}
# Get the names for the output data set
n1 <- length(eff$snp.levels)
n2 <- length(eff$var.levels)
if (flag) n2 <- n2 + 1
names <- matrix(data=" ", nrow=n1, ncol=n2)
if (!flag) {
for (i in 1:n1) {
names[i,] <- paste(eff$var1, eff$snp.levels[i], "_",
eff$var2, "_", eff$var.levels, sep="")
}
} else {
temp <- c(baseline, eff$var2)
for (i in 1:n1) {
names[i,] <- paste(eff$var1, eff$snp.levels[i], "_",
temp, sep="")
}
}
eff$out.names <- names
# Get the list names
temp <- character(length(eff$type))
i <- 1
for (type in eff$type) {
if (type == 1) {
temp[i] <- ".JointEffects"
} else {
temp[i] <- ".StratEffects"
}
i <- i + 1
}
eff$list.names <- temp
eff
} # END: setupEffects
# Check the input lists
snp.list <- check.snp.list(snp.list)
pheno.list <- default.list(pheno.list,
c("response.var"), list("ERROR"), error=c(1))
# Rename main.vars to X.vars, int.vars to V.vars
pheno.list$X.vars <- pheno.list[["main.vars", exact=TRUE]]
pheno.list$V.vars <- pheno.list[["int.vars", exact=TRUE]]
# Determine if fomulas were passed in
temp <- pheno.list[["X.vars", exact=TRUE]]
if (!is.null(temp)) {
main.form <- ("formula" %in% class(temp))
} else {
main.form <- FALSE
}
temp <- pheno.list[["V.vars", exact=TRUE]]
if (!is.null(temp)) {
int.form <- ("formula" %in% class(temp))
} else {
int.form <- FALSE
}
temp <- pheno.list[["strata.var", exact=TRUE]]
if (!is.null(temp)) {
s.form <- ("formula" %in% class(temp))
} else {
s.form <- FALSE
}
formFlag <- main.form + int.form + s.form
temp.list <- op[["temp.list", exact=TRUE]]
temp.list <- check.temp.list(temp.list)
temp <- c("BFGS", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "IRLS")
op <- default.list(op,
c("id", "print", "optimizer", "snpName", "tests.CML",
"tests.EB", "genetic.model", "geno.counts", "subject.counts",
"test.omnibus", "test.main", "test.inter", "tests.UML",
"allele.cc"),
list(1, 0, "BFGS", "SNP_", 1, 1, 0, 0, 0, 0, 0, 0, 1, 1),
checkList=list(NA, 0:1, temp, NA, 0:1, 0:1, 0:3,
0:1, 0:1, 0:1, 0:1, 0:1, NA, 1:2))
op <- checkOptions(op)
op$tests.UML <- check.tests.uml(getListName(op, "tests.UML"))
# Check cc.var
temp <- pheno.list[["cc.var", exact=TRUE]]
if (is.null(temp)) {
op$subject.counts <- 0
if (op$allele.cc == 2) pheno.list$cc.var <- pheno.list$response.var
}
# Check the out.est list
out.est <- op[["out.est", exact=TRUE]]
if (!is.null(out.est)) {
out.est <- check.out.est(out.est)
out.est.flag <- 1
# Set flags for output
out.est.beta <- ("beta" %in% out.est$what)
out.est.se <- ("se" %in% out.est$what)
out.est.test <- ("test" %in% out.est$what)
out.est.pval <- ("pvalue" %in% out.est$what)
out.est.corr <- !is.null(getListName(out.est, "corr.parms"))
} else {
out.est.flag <- 0
}
# The genetic model is taken care of in logistic.SNP
snp.list$genetic.model <- NULL
# Flag for streamed input
stream <- snp.list[["stream", exact=TRUE]]
# All the variables must be given by variable name (not column number)
if ((is.numeric(pheno.list$response.var)) ||
(is.numeric(pheno.list$strata.var)) ||
(is.numeric(pheno.list$factor.vars)) ||
(is.numeric(pheno.list$id.var)) ) {
stop("ERROR: variables must be specified by name, not column number")
}
if ((!main.form) && (is.numeric(pheno.list$X.vars)))
stop("ERROR: variables must be specified by name, not column number")
if ((!int.form) && (is.numeric(pheno.list$V.vars)))
stop("ERROR: variables must be specified by name, not column number")
# Keep only the variables we need
if (!formFlag) {
temp <- c(pheno.list$response.var, pheno.list$strata.var,
pheno.list$X.vars, pheno.list$V.vars, pheno.list$id.var,
pheno.list$cc.var)
pheno.list$keep.vars <- unique(temp)
pheno.list$remove.vars <- NULL
} else {
# Do not remove variables
pheno.list$keep.vars <- NULL
}
pheno.list$remove.miss <- 1
pheno.list$make.dummy <- 0
# Get the data vector of snps
tlist <- list(include.row1=0, include.snps=0, return.type=1, MAF=1,
missing=1, snpNames=1, orderByPheno=1, return.pheno=1)
if (stream) {
tlist$file.index <- 1
temp <- paste("TMP_", temp.list$id, sep="")
tempfile <- getTempfile(temp.list$dir, temp)
tlist$outfile <- tempfile
}
temp <- try(getData.1(snp.list, pheno.list, temp.list, op=tlist),
silent=TRUE)
if (class(temp) == "try-error") {
print(temp)
stop("ERROR loading data")
}
snpData <- temp$data
missing <- temp$missing
snpNames <- temp$snpNames
delimiter <- getDelimiter(snp.list, output=1)
design.X <- NULL
design.V <- NULL
nsnps <- length(snpData)
maf.vec <- temp[["MAF", exact=TRUE]]
maf.flag <- !is.null(maf.vec)
controls <- temp[["controls", exact=TRUE]]
alleles <- temp[["alleles", exact=TRUE]]
allFlag <- !is.null(alleles)
ProbG1 <- temp[["ProbG1", exact=TRUE]]
ProbG1.flag <- !is.null(ProbG1)
# Get the phenotype data
phenoData.list <- temp$phenoData.list
phenoData0 <- phenoData.list$data
if (stream) {
snpFid <- temp$fid
in.miss <- getListName(snp.list, "in.miss")
heter.codes <- getListName(snp.list, "heter.codes")
delete <- temp.list$delete
# Update the snp list
snp.list <- update.snp.list(snp.list, where=1)
snpNames <- getListName(snp.list, "snpNames")
snpFlag <- 1 - is.null(snpNames)
# Get the required objects
temp <- setUp.stream(temp, snp.list, tempfile, delete, controls=controls)
subjFlag <- temp$subjFlag
subj.ids <- temp$subj.ids
total.nsubjects <- temp$total.nsubjects
subj.order2 <- temp$subj.order2
start.stream <- temp$start.stream
stop.stream <- temp$stop.stream
snp <- temp$snp
control.ids <- getListName(temp, "control.ids")
codes <- getInheritanceVec(snp.list$genetic.model,
recode=snp.list$recode)
}
rm(temp, tlist, snp.list, phenoData.list, controls)
temp <- gc(verbose=FALSE)
# Get the response variable
response0 <- as.numeric(phenoData0[, pheno.list$response.var])
if (!any(response0 %in% 0:1)) stop("ERROR: response must be 0-1")
# Get the number of observations
nobs <- length(response0)
# Print out number of observations that will be used
temp <- paste("For the analysis, ", nobs, " observations will be used.", sep="")
print(temp)
print(table(response0))
# Determine all factor variables
facVars <- pheno.list[["factor.vars", exact=TRUE]]
Xflag <- !is.null(pheno.list[["X.vars", exact=TRUE]])
Vflag <- !is.null(pheno.list[["V.vars", exact=TRUE]])
Sflag <- !is.null(pheno.list[["strata.var", exact=TRUE]])
vars <- pheno.list$response.var
if ((!s.form) && (Sflag)) vars <- c(vars, pheno.list$strata.var)
if ((!main.form) && (Xflag)) vars <- c(vars, pheno.list$X.vars)
if ((!int.form) && (Vflag)) vars <- c(vars, pheno.list$V.vars)
vars <- unique(vars)
for (v in vars) {
if ((is.factor(phenoData0[, v])) || (is.character(phenoData0[, v]))) {
facVars <- c(facVars, v)
}
}
facVars <- unique(facVars)
facFlag <- !is.null(facVars)
X.newVars <- NULL
V.newVars <- NULL
# Get the stratafication matrix
temp <- sMatrix.logistic(phenoData0, pheno.list$strata.var, facVars)
design.S0 <- temp$design.S0
op$s.1catVar <- temp$s.1catVar
nStrata <- ncol(design.S0)
# Get the X design matrix
if (Xflag) {
temp <- logistic.dsgnMat(phenoData0, pheno.list$X.vars,
facVars)
design.X0 <- temp$designMatrix
X.newVars <- temp$newVars
} else {
design.X0 <- NULL
}
# Get the V design matrix
if (Vflag) {
temp <- logistic.dsgnMat(phenoData0, pheno.list$V.vars,
facVars)
design.V0 <- temp$designMatrix
V.newVars <- temp$newVars
} else {
design.V0 <- NULL
op$test.inter <- 0
}
# Get the genetic model
genetic.model <- op$genetic.model
# Check for effects
op$effects <- setupEffects(getListName(op, "effects"), facVars,
pheno.list$X.vars, pheno.list$V.vars)
effects <- getListName(op, "effects")
effectsFlag <- !is.null(effects)
# Get a logical vector for obtaining the cc counts
if (op$subject.counts) {
cntrl <- (pheno.list$cc.var == 0)
subj.cnts0 <- table(cntrl)
}
rm(pheno.list, phenoData0)
temp <- gc(verbose=FALSE)
# Set the output directory
if (!is.null(op$out.dir)) {
out.dir <- checkForSep(op$out.dir)
if (nsnps != length(unique(snpNames))) stop("SNP names are not unique")
} else {
out.dir <- NULL
}
# Determine if tests will be computed
omniFlag <- op$test.omnibus
mainFlag <- op$test.main
interFlag <- op$test.inter
test2Flag <- (!is.null(getListName(op, "tests")))
tests.UML <- (!is.null(getListName(op, "tests.UML")))
if (tests.UML) tests.UML.vec <- op$tests.UML
tests.CML <- op$tests.CML
tests.EB <- op$tests.EB
tests <- omniFlag + mainFlag + test2Flag + tests.UML +
tests.CML + tests.EB + interFlag
geno.counts <- op$geno.counts
subj.counts <- op$subject.counts
# Check out.est
if (out.est.flag) {
if ((out.est$parms == 2) && (!Vflag)) out.est$parms <- 1
}
# Get the variable names that will be used
log.vnames <- logistic.vnames(design.X0, design.V0, nStrata,
snpName=op$snpName, out.est=out.est,
genetic.model=genetic.model)
est.p <- log.vnames[["est.p", exact=TRUE]]
est.parms <- log.vnames[["est.parms", exact=TRUE]]
# Get the variable names for the tests
if (mainFlag) main.vars <- getMain.vars()
if (omniFlag) omni.vars <- getOmni.vars()
if (interFlag) inter.vars <- getInter.vars()
if (test2Flag) {
# test2.vars is a list
test2.vars <- getTest2.vars()
op$test2.vars <- NULL
n.tests <- length(test2.vars)
}
rm(facFlag, facVars, X.newVars, V.newVars)
# Open the output file
fileFlag <- 1 - is.null(op$out.file)
if (fileFlag) {
temp <- initOutfile()
fid <- temp$fid
outVec <- temp$outVec
est.se <- temp$est.se
est.test <- temp$est.test
est.pval <- temp$est.pval
est.corr <- temp$est.corr
}
rm(log.vnames)
# Base model
temp <- getListName(op, "base.outfile")
if (!is.null(temp)) {
i <- callGLM(response0, X.main=design.X0, X.int=design.V0, int.vec=NULL,
family=binomial(), prefix="SNP_", retX=TRUE,
retY=TRUE, inc.int.vec=0, int.vec.base=0)
save(i, file=temp)
print(summary(i))
}
op$errorCheck <- 0
ProbG1.vec <- NULL
if (ProbG1.flag) op$imputed <- 1
# Loop over each SNP
i <- 0
while (1) {
if (fileFlag) outVec[] <- NA
i <- i + 1
if (!stream) {
if (i > nsnps) break
snp <- as.numeric(getVecFromStr(snpData[i], delimiter=delimiter))
snp.miss <- missing[i]
snp.name <- snpNames[i]
if (allFlag) {
majMin <- alleles[i]
} else {
majMin <- " "
}
# Get the MAF
if (maf.flag) {
maf <- maf.vec[i]
} else {
maf <- getMAF(snp)
}
# For imputed data
if (ProbG1.flag) ProbG1.vec <- as.numeric(getVecFromStr(ProbG1[i], delimiter=delimiter))
} else {
if (i > 1) {
start.stream <- start.stream + 1
snp <- getNextObs(i, snpFid, snpFlag, snpNames, tempfile, delete,
start.stream, stop.stream, delimiter)
if (is.null(snp)) break
}
snp.name <- snp[1]
if (snpFlag) snpNames <- update.snpNames(snpNames, snp.name)
# Get the correct order and recode
temp <- orderSNP(snp[-1], snp.name, subj.order2, total.nsubjects,
in.miss, heter.codes, subjFlag=subjFlag, subj.ids=subj.ids,
out.genotypes=codes, control.ids=control.ids)
snp <- as.numeric(temp$SNP)
maf <- temp$MAF
snp.miss <- any(is.na(snp))
if (allFlag) {
majMin <- temp$alleles
} else {
majMin <- " "
}
}
# Remove missing values
if (snp.miss) {
temp <- as.logical(!is.na(snp))
if (subj.counts) subj.notMiss <- temp
response <- response0[temp]
snp <- snp[temp]
design.S <- removeOrKeepRows(design.S0, temp, which=1)
if (Xflag) design.X <- removeOrKeepRows(design.X0, temp, which=1)
if (Vflag) design.V <- removeOrKeepRows(design.V0, temp, which=1)
} else {
response <- response0
design.S <- design.S0
design.X <- design.X0
design.V <- design.V0
subj.notMiss <- NULL
}
# Get the genotype counts
genoCounts <- getGenoCounts(snp)
n.genos <- sum(genoCounts > 0)
# Make sure that there is more than 1 genotype
if (n.genos < 2) {
# Do not call snp.logistic
ret <- list()
} else {
# Set options
if (n.genos == 2) {
op$geno.binary <- 1
op$genetic.model <- 0
} else {
op$geno.binary <- 0
op$genetic.model <- genetic.model
}
# Fit the model
ret <- try(snp.main(response, snp, X.main=design.X, X.int=design.V,
X.strata=design.S, ProbG1=ProbG1.vec, op=op), silent=TRUE)
if ("try-error" %in% class(ret)) ret <- list()
}
# Get the list index for test variable names
vListIndex <- getVListIndex(n.genos)
# Add the snp name
ret$snp <- snp.name
# Compute tests
if (tests) {
if (tests.UML) ret <- addToList(ret, "UML")
if (tests.CML) ret <- addToList(ret, "CML")
if (tests.EB) ret <- addToList(ret, "EB")
}
# Effects
if (effectsFlag) ret <- calcEffects(effects, ret)
# Print the return list
if (op$print) print(ret)
# Save the output
if (!is.null(out.dir)) {
temp <- paste(out.dir, "out_", snp.name, ".rda", sep="")
save(ret, file=temp)
}
if (fileFlag) temp <- outputRow()
} # END: while(1)
if (fileFlag) close(fid)
ret
} # END: snp.scan.logistic
# Function to permform a SNP by environment interaction analysis for 1 SNP.
snp.logistic <- function(data, response.var, snp.var, main.vars=NULL,
int.vars=NULL, strata.var=NULL, op=NULL) {
# INPUT:
# data Data frame containing all the data.
# No default
# response.var Name of the binary response variable coded
# as 0 (controls) and 1 (cases)
# No default.
# snp.var Name of the genotype variable coded as 0, 1, 2 (or 0, 1)
# No default.
# main.vars Character vector of the main effects variables or
# a formula.
# The default is NULL
# int.vars Names of all covariates of interest that will
# interact with the SNP variable or a formula.
# The default is NULL
# strata.var Name of the stratification variable(s) or a formula.
# Any character variable or factor will be considered
# categorical. Any numeric variable will be considered
# continuous. An intercept column will be automatically
# included.
# The default is NULL (1 stratum)
# ProbG1.var Variable for Prob(G = 1) or NULL. Not needed if snp.var
# is of length 3.
# The default is NULL.
#####################################################################
# op List with the following names.
# genetic.model 0-2
# 0: trend
# 1: dominant
# 2: recessive
# 3: general
# The default is 0.
# This option has no effect if the snp is binary.
# reltol Stopping tolerance
# The default is 1e-8
# maxiter Maximum number of iterations
# The default is 100
# optimizer One of :"Nelder-Mead", "BFGS", "CG", "L-BFGS-B",
# "SANN", "IRLS"
# The default is "BFGS"
# snpName Name to be used in the SNP variable and interaction
# variables.
# The default is "SNP_".
# debug 0 or 1 to show the IRLS iterations.
# The default is 0
# fit.null 0 or 1 to fit a NULL model which excludes the snp,
# interactions with the snp, and the interacting covariates
# This option takes precedence over zero.vars.
# The default is 0.
# zero.vars Character vector of variable names to fix or a list
# with names "snp.var", "main.vars", int.vars", and
# "strata.var".
# If zero.vars is a list, then the names is the list can
# either be formulas or character vectors.
# The default is 0.
######################################################################
# Check for errors
if (length(response.var) != 1) stop("response.var must be a single variable")
if (!is.data.frame(data)) stop("data must be a data frame")
ProbG1.var <- NULL
n.snp.var <- length(snp.var)
#if (!(n.snp.var %in% c(1, 3))) {
# stop("snp.var must be a single variable or 3 variables for imputed genotypes")
#}
if (n.snp.var != 1) stop("snp.var must be a single variable")
# Check variable names
vlist <- list(response.var=response.var, snp.var=snp.var, main.vars=main.vars,
int.vars=int.vars, strata.var=strata.var, ProbG1.var=ProbG1.var)
vars <- getAllVars(vlist, names=names(vlist))
temp <- !(vars %in% colnames(data))
if (any(temp)) {
print(vars[temp])
stop("The above variables were not found in the input data")
}
# Check if snp.var is in main.vars or int.vars
if (any(snp.var %in% getAllVars(main.vars))) stop("ERROR: main.vars must not contain snp.var")
if (any(snp.var %in% getAllVars(int.vars))) stop("ERROR: int.vars must not contain snp.var")
main.call <- main.vars
int.call <- int.vars
strata.call <- strata.var
op <- default.list(op, c("snpName", "fit.null", "imputed", "genetic.model"),
list("SNP_", 0, 0, 0))
op$imputed <- 0
if ((!is.numeric(snp.var)) && (n.snp.var == 1)) op$snpName <- snp.var
zeroFlag <- 0
zero.vars <- NULL
#if ((n.snp.var > 1) || (!is.null(ProbG1.var))) {
# op$imputed <- 1
#}
if (n.snp.var == 1) {
snp <- as.numeric(unfactor(data[, snp.var]))
snp <- snp[is.finite(snp)]
#if (!all(snp %in% 0:2)) op$imputed <- 1
if (!all(snp %in% 0:2)) stop("snp.var must be coded 0-1-2")
}
imputed <- op$imputed
genetic.model <- op$genetic.model
if (!(genetic.model %in% 0:3)) stop("op$genetic.model must be 0-3")
#if ((imputed) && (genetic.model == 3)) stop("op$genetic.model must be 0-2 for imputed genotypes")
#if ((imputed) && (n.snp.var == 1) && (genetic.model != 0)) {
# stop("op$genetic.model must be 0 for imputed genotypes and length(snp.var) = 1")
#}
#if ((imputed) && (n.snp.var == 1) && (is.null(ProbG1.var))) {
# stop("ProbG1.var must be specified for imputed genotypes and length(snp.var) = 1")
#}
#if ((imputed) && (length(ProbG1.var) > 1)) stop("ProbG1.var must be NULL or of length 1")
main.form <- ("formula" %in% class(main.vars))
int.form <- ("formula" %in% class(int.vars))
s.form <- ("formula" %in% class(strata.var))
ProbG1 <- NULL
# For impute snp with 3 vars, if all 0, then missing
#if (n.snp.var > 1) {
# snp1 <- as.numeric(unfactor(data[, snp.var[1]]))
# snp2 <- as.numeric(unfactor(data[, snp.var[2]]))
# snp3 <- as.numeric(unfactor(data[, snp.var[3]]))
# temp <- (snp1 %in% 0) & (snp2 %in% 0) & (snp3 %in% 0)
# if (any(temp)) data[temp, snp.var] <- NA
#}
# Remove missing values
temp <- getFormulas(vlist)
miss <- c(NA, NaN, Inf, -Inf)
if (length(temp)) data <- applyFormulas(data, temp, remove=miss)
data <- removeMiss.vars(data, vars=vars, miss=miss)
rm(vlist, vars)
temp <- gc(verbose=FALSE)
# Get the response variable
D <- unfactor(data[, response.var])
nobs <- length(D)
# Get Prob(G=1)
#if ((imputed) && (!is.null(ProbG1.var))) ProbG1 <- as.numeric(unfactor(data[, ProbG1.var]))
# Get the snp variable(s)
if (n.snp.var == 1) {
snp <- as.numeric(unfactor(data[, snp.var]))
} else {
# snp1 <- as.numeric(unfactor(data[, snp.var[1]]))
# snp2 <- as.numeric(unfactor(data[, snp.var[2]]))
# snp3 <- as.numeric(unfactor(data[, snp.var[3]]))
# if (genetic.model == 0) {
# snp <- snp2 + 2*snp3
# } else if (genetic.model == 1) {
# snp <- snp2 + snp3
# } else if (genetic.model == 2) {
# snp <- snp3
# }
# Add it ot the data frame
# data[, op$snpName] <- snp
# Save Prob(G = 1) if needed
# if (is.null(ProbG1)) ProbG1 <- snp2
# rm(snp1, snp2, snp3)
# gc()
}
#if (imputed) op$genetic.model <- 0
facVars <- NULL
sflag <- !is.null(strata.var)
# Check for constant strata variable
if ((sflag) && (!s.form) && (length(strata.var) == 1)) {
temp <- unique(makeVector(data[, strata.var]))
if (length(temp) == 1) sflag <- FALSE
}
# Determine if any strata vars are character
if ((sflag) && (!s.form)) {
for (v in strata.var) {
if (is.character(data[, v])) data[, v] <- factor(data[, v])
}
}
# Get the variables that are factors
for (temp in colnames(data)) {
if (is.factor(data[, temp])) facVars <- c(facVars, temp)
}
# Get the V design matrix
design.V0 <- logistic.dsgnMat(data, int.vars, facVars)$designMatrix
int.vars <- colnames(design.V0)
# For fitting a null model
if (op$fit.null) {
# Remove interaction vars, snp
temp <- getAllVars(int.call)
temp <- unique(temp, int.vars)
zero.vars <- temp
int.vars <- NULL
design.V0 <- NULL
zeroFlag <- 1
op$fixParms <- list(parms=snp.var, values=0)
op$zero.vars <- NULL
}
# Get the stratafication matrix
temp <- sMatrix.logistic(data, strata.var, facVars)
design.S0 <- temp$design.S0
op$s.1catVar <- temp$s.1catVar
# Get the X design matrix
design.X0 <- logistic.dsgnMat(data, main.vars, facVars, remove.vars=zero.vars)$designMatrix
# For the zero.vars option
zero.vars <- op[["zero.vars", exact=TRUE]]
if (!is.null(zero.vars)) {
temp <- list(main.vars=design.X0, int.vars=design.V0, strata.var=design.S0)
temp <- apply_zero.vars(zero.vars, temp, snp.var, facVars, data)
op$fixParms <- temp[["fixParms", exact=TRUE]]
temp <- temp$mat.list
design.V0 <- temp[["int.vars", exact=TRUE]]
int.vars <- colnames(design.V0)
design.X0 <- temp[["main.vars", exact=TRUE]]
design.S0 <- temp[["strata.var", exact=TRUE]]
}
main.vars <- colnames(design.X0)
strata.var <- colnames(design.S0)
# Check for non-dummy continuous variables ant print a warning
wflag <- 0
if (!all(design.X0 %in% 0:1)) wflag <- 1
if ((!wflag) && (!all(design.V0 %in% 0:1))) wflag <- 1
if ((!wflag) && (!all(design.S0 %in% 0:1))) wflag <- 1
if (wflag) warning("For stability of the algorithm, continuous variables should be scaled")
# Call the core function
ret <- snp.main(D, snp, X.main=design.X0, X.int=design.V0,
X.strata=design.S0, ProbG1=ProbG1, op=op)
# Add model info
model <- list(data=data, response.var=response.var, snp.var=op$snpName,
main.vars=main.vars, int.vars=int.vars,
strata.var=strata.var, factors=facVars,
snpName=op$snpName, main.call=main.call,
int.call=int.call, strata.call=strata.call)
ret$model.info <- model
ret
} # END: snp.logistic
# Function to permform a SNP by environment interaction analysis for 1 SNP.
snp.main <- function(D, snp, X.main=NULL, X.int=NULL,
X.strata=NULL, ProbG1=NULL, op=NULL) {
# INPUT:
# D Binary response vector coded as 0 (controls) and 1 (cases)
# No default.
# snp Vector of genotypes or 3 column matrix for imputed snp
# No default.
# X.main Design matrix for all covariates of interest, excluding
# the SNP variable. This matrix should NOT include an
# intercept column.
# The default is NULL
# X.int Design matrix for all covariates of interest that will
# interact with the SNP variable.
# This matrix should NOT include an intercept column.
# The default is NULL
# X.strata NULL or a design matrix for the stratification.
# The default is NULL (1 stratum)
# ProbG1 Vector of Prob(G=1) for imputed SNPs
# The default is NULL
#####################################################################
# op List with the following names.
# genetic.model 0-3
# 0: trend
# 1: dominant
# 2: recessive
# 3: general
# The default is 0.
# This option has no effect if the snp is binary.
# reltol Stopping tolerance
# The default is 1e-8
# maxiter Maximum number of iterations
# The default is 100
# optimizer One of :"Nelder-Mead", "BFGS", "CG", "L-BFGS-B",
# "SANN"
# The default is "BFGS"
# snpName Name to be used in the SNP variable and interaction
# variables.
# The default is "SNP_".
# debug 0 or 1 to show the IRLS iterations.
# The default is 0
# errorCheck 0 or 1 to check for errors with the input arguments
# The default is 1
# use.C.code 0 or 1 to call the C code for the BFGS optimizer
# The default is 1
# fit.null 0 or 1 for fitting a null model. All the necessary
# variables should already be removed from the design
# matrices.
# The default is 0
# num.deriv 0 or 1 to use numerical derivatives.
# Numerical derivatives will be used if genetic.model > 0
# The default is 0
# imputed 0 or 1 for imputed SNPs
# s.1catVar 0 or 1 for 1 categorical strata var
# The default is 0
###################################################################
# fixParms List with names "parms" and "value".
# The default is NULL
# parms Vector of parm names to fix
# No default
# values Numeric vector of fixed values
# No default
######################################################################
# RETURN:
# A list with sublists named UML (standard logistic regression),
# CML (conditional ML), and EB (empirical bayes).
# Each sublist contains the parameter estimates and covariance
# matrix. The lists UML and CML also contain the log-likelihood.
# Local function to return the initial estimates
getInit <- function() {
# Define the formula
ff <- "D ~ 1"
if (nx) ff <- paste(ff, " + X.main ", sep="")
if (!fit.null) {
if (!zeroSNP) ff <- paste(ff, "+ fsnp ", sep="")
if (nv) {
if (!gmodel3) {
ff <- paste(ff, "+ fsnp:X.int", sep="")
} else {
ff <- paste(ff, "+ fsnp[,1]:X.int + fsnp[,2]:X.int", sep="")
}
}
}
ff <- as.formula(ff)
# Get variable names for Z
v <- logistic.vnames(X.main, X.int, nStrata, snpName=op$snpName,
genetic.model=genetic.model, fit.null=fit.null, zeroSNP=zeroSNP)
# Call after snp was multiplied by V
fit <- glm(ff, family=binomial(), model=FALSE, x=FALSE, y=FALSE)
# Check the convergence
if (!fit$converged) return(NULL)
UML.parms <- fit$coefficients
loglike <- getLoglike.glm(fit)
cov <- summary(fit)$cov.scaled
cnames <- v$UML.names
#all <- v$all.names
names(UML.parms) <- cnames
parms <- UML.parms
if (zeroSNP) {
# Add the snp back in
pos <- nx + 1
len <- length(parms)
temp <- c(parms[1:pos], 0)
temp2 <- c(cnames[1:pos], v$snp)
if (nv) {
temp <- c(temp, parms[(pos+1):len])
temp2 <- c(temp2, cnames[(pos+1):len])
}
parms <- temp
names(parms) <- temp2
}
# Check for NAs
naFlag <- 0
naXPos <- NULL
naVPos <- NULL
temp <- is.na(UML.parms)
cnames <- cnames[!temp]
colnames(cov) <- cnames
rownames(cov) <- cnames
if (any(temp)) {
naFlag <- 1
# Stop if snp was not estimated
if (temp[nx+2]) stop("ERROR: SNP was not estimated for the UML method")
# Get the NA positions in the design matrices
naPos <- (1:length(UML.parms))[temp]
temp2 <- (naPos <= nx + 1)
if (any(temp2)) naXPos <- naPos[temp2] - 1
temp2 <- !temp2
if (any(temp2)) naVPos <- naPos[temp2] - nx - 2
UML.parms <- UML.parms[!temp]
parms <- parms[!is.na(parms)]
}
# Remove the special character "`" that glm will sometimes
# put in the variable names
parms <- changeStr.names(parms, "`", replace="")
cov <- changeStr.names(cov, "`", replace="")
alpha <- parms[1]
beta <- parms[-1]
# Vector for xi parms. Only intercept will be non zero
xi <- rep(0, times=nStrata)
total <- 2*length(D)
if (op$s.1catVar) {
mini <- NULL
maxi <- NULL
for (i in 1:ncol(X.strata)) {
temp <- X.strata[, i] == 1
freq <- sum(snp[temp])/total
if (freq == 0) {
mini <- c(mini, i)
} else if (freq == 1) {
maxi <- c(maxi, i)
} else {
xi[i] <- log(freq/(1-freq))
}
if (!is.finite(xi[i])) xi[i] <- 0
if (!is.null(mini)) xi[mini] <- min(xi, -5)
if (!is.null(maxi)) xi[maxi] <- max(xi, 5)
}
} else {
freq <- sum(snp)/total
if (freq == 0) {
xi[1] <- -1
} else if (freq == 1) {
xi[1] <- 1
} else {
xi[1] <- log(freq/(1-freq))
}
if (!is.finite(xi[1])) xi[1] <- 0
}
eta <- c(parms, xi)
names(eta) <- c(names(parms), v$strata)
list(eta=eta, alpha=alpha, beta=beta, xi=xi, fit=fit, parms=parms,
loglike=loglike, cov=cov, fitted.values=fit$fitted.values, UML.parms=UML.parms,
naFlag=naFlag, naXPos=naXPos, naVPos=naVPos)
} # END: getInit
# Function to compute Pdg.xs = P(D=d, G=g | X, S)
Pdg.xs <- function(ret, alpha, beta, xi) {
# Get the xi parameters for each observation
if (nStrata == 1) {
temp.xi <- xi
} else {
dim(xi) <- c(nStrata, 1)
temp.xi <- X.strata %*% xi
}
# Make sure that beta is a column vector
dim(beta) <- c(nbeta, 1)
# Get theta for each d, g combination
ret[, d0g0.col] <- 0
ret[, d0g1.col] <- log2 + temp.xi
ret[, d0g2.col] <- 2*temp.xi
ret[, d1g0.col] <- alpha + (Z0 %*% beta)
ret[, d1g1.col] <- alpha + (Z1 %*% beta) + log2 + temp.xi
ret[, d1g2.col] <- alpha + (Z2 %*% beta) + 2*temp.xi
ret <- exp(ret)
# For a binary snp that the user input
if (geno.binary) ret[, g2.col] <- 0
# Compute the sum over (d,g)
sum <- rowSums(ret)
# Divide by sum
ret <- matrixDivideVec(ret, sum)
list(Pdg.matrix=ret, Pdg.rowSums=sum)
} # END: Pdg.xs
# Function to compute P(D=1 | X,S)
Pd1.xs <- function(pmat) {
rowSums(pmat[, d1.col])
} # END: Pd1.xs
# Function to compute E(DG | X,S)
Edg.xs <- function(pmat) {
if (gmodel3) {
return(pmat[, c(d1g1.col, d1g2.col)])
} else if (geno.binary) {
return(pmat[, d1g1.col])
} else {
return(pmat[, d1g1.col] + 2*pmat[, d1g2.col])
}
} # END: Edg.xs
# Function to compute E(G | X,S)
Eg.xs <- function(pmat) {
temp <- pmat
temp[, g2.col] <- 2*temp[, g2.col]
temp <- temp[, c(g1.col, g2.col)]
rowSums(temp)
} # END: Eg.xs
# Function to compute E(DG^2 | X,S) = E(D^2G^2 | X,S)
Edgg.xs <- function(pmat) {
if (geno.binary) {
return(pmat[, d1g1.col])
} else {
return(pmat[, d1g1.col] + 4*pmat[, d1g2.col])
}
} # END: Edgg.xs
# Function to compute E(G^2 | X,S)
Egg.xs <- function(pmat) {
temp <- pmat
temp[, g2.col] <- 4*temp[, g2.col]
temp <- temp[, c(g1.col, g2.col)]
rowSums(temp)
} # END: Egg.xs
# Function to return a logical matrix to make the computation of
# the log-likelihood easier.
getLoglike.mat <- function() {
ret <- matrix(data=FALSE, nrow=n, ncol=nlevels)
if (gmodel3) {
col <- 3*D + 1 + snp
} else {
col <- 3*D + 1 + fsnp
}
for (i in 1:n) ret[i, col[i]] <- TRUE
ret
} # END: getLoglike.mat
# Function to compute the log-likelihood (or testing)
loc.getLoglike <- function(eta) {
alpha <- eta[alpha.row]
beta <- eta[beta.row]
xi <- eta[xi.row]
if (fixFlag) eta <- fixGetEta(eta)
Pdg <- Pdg.xs(Pdg, eta[alpha.row], eta[beta.row], eta[xi.row])
if (!imputed) {
Pdg <- Pdg$Pdg.matrix
ret <- sum(log(Pdg[loglike.mat]))
} else {
# Get the xi parameters for each observation
if (nStrata == 1) {
temp.xi <- xi
} else {
dim(xi) <- c(nStrata, 1)
temp.xi <- X.strata %*% xi
}
# Make sure that beta is a column vector
dim(beta) <- c(nbeta, 1)
vec <- D*(alpha + Z.imp %*% beta) + fsnp*temp.xi + log2*ProbG1
vec <- exp(vec)/Pdg$Pdg.rowSums
ret <- sum(log(vec))
}
ret
} # END: loc.getLoglike
# Function to compute a Z matrix
getZ <- function(gvalue) {
temp.V <- NULL
if (gmodel3) {
temp.snp <- matrix(data=0, nrow=n, ncol=2)
if (gvalue) temp.snp[, gvalue] <- 1
if (nv) {
temp.V <- cbind(matrixMultVec(X.int, temp.snp[, 1]),
matrixMultVec(X.int, temp.snp[, 2]))
}
} else {
temp.snp <- rep(gvalue, times=n)
if (nv) temp.V <- matrixMultVec(X.int, temp.snp)
}
ret <- as.matrix(cbind(X.main, temp.snp, temp.V))
ret
} # END: getZ
# Function to compute Z matrix for imputed snp
getZ.imp <- function(genovec) {
temp.V <- NULL
if (nv) temp.V <- matrixMultVec(X.int, genovec)
ret <- as.matrix(cbind(X.main, genovec, temp.V))
ret
} # END: getZ.imp
# Function to compute W(Y - mu)
getWtYmu <- function(eta=NULL, which=1) {
# eta Vector of parms
# Only needed for which = 1
# which 0 or 1, 1 is for the optimizer
# The default is 1
# Get the matrix of probabilities
if (fixFlag) eta <- fixGetEta(eta)
pmat <- Pdg.xs(Pdg, eta[alpha.row], eta[beta.row], eta[xi.row])$Pdg.matrix
temp1 <- D - Pd1.xs(pmat)
temp2 <- Dsnp - Edg.xs(pmat)
temp3 <- snp - Eg.xs(pmat)
dim(temp1) <- c(1, n)
temp1 <- temp1 %*% X.main
if (gmodel3) {
dim(temp2) <- c(2, n)
} else {
dim(temp2) <- c(1, n)
}
temp2 <- temp2 %*% X.int
dim(temp3) <- c(1, n)
temp3 <- temp3 %*% X.strata
temp <- c(temp1, temp2, temp3)
dim(temp) <- c(nparms, 1)
if (fixFlag) temp <- temp[fixMap]
temp
} # END: getWtYmu
# Function to call the optimizer
callOptim <- function() {
# Set the control list
control <- list(fnscale=-1, maxit=op$maxiter, reltol=op$reltol)
# Call the optimizer
if (op$num.deriv) {
ret <- optim(eta0, loc.getLoglike, method=op$optimizer,
control=control, hessian=TRUE)
} else {
ret <- optim(eta0, loc.getLoglike, gr=getWtYmu, method=op$optimizer,
control=control, hessian=TRUE)
}
# Determine if it converged
conv <- (ret$convergence == 0)
# Get the covariance matrix and score
if (conv) {
conv <- 1
cov <- chol(-ret$hessian)
cov <- chol2inv(cov)
cnames <- names(eta0)
colnames(cov) <- cnames
rownames(cov) <- cnames
} else {
cov <- NA
conv <- 0
}
# Return list
list(parms=ret$par, cov=cov, converged=conv, loglike=ret$value)
} # END: callOptim
# Function to compute empirical bayes estimates
getEB <- function() {
#if (fixFlag) return(NULL)
ids <- c(alpha.row, beta.row)
if (zeroSNP) ids <- ids[-length(ids)]
parm1 <- ret$UML$parms
parm2 <- ret$CML$parms[ids]
cov1 <- ret$UML$cov
vnames <- names(parm1)
dim(parm1) <- NULL
dim(parm2) <- NULL
psi2 <- parm1 - parm2
psi2 <- psi2*psi2
phi <- diag(cov1)
denom <- psi2 + phi
# Compute the new estimates
parms <- parm1*psi2/denom + parm2*phi/denom
temp <- (phi*(phi-psi2))/(denom*denom)
if (length(temp) > 1) {
cmat <- cbind(diag(1-temp), diag(temp))
} else {
cmat <- matrix(c(1-temp, temp), nrow=1)
}
rm(psi2, phi, denom, parm1, parm2)
temp <- gc(verbose=FALSE)
# Set up the Z matrix
temp <- NULL
if (nv) {
temp <- removeOrKeepCols(X.int, 1, which=-1)
if (gmodel3) {
temp <- cbind(matrixMultVec(temp, fsnp[, 1]),
matrixMultVec(temp, fsnp[, 2]))
} else {
temp <- matrixMultVec(temp, fsnp)
}
}
# Define the Z matrix. X has an intercept column as the first column
if (!zeroSNP) {
Z <- cbind(X.main, fsnp, temp)
} else {
Z <- cbind(X.main, temp)
}
# Compute the scores for UML
temp <- D - fitted.values
score1 <- matrixMultVec(Z, temp)
# Get the matrix of probabilities P(D=d, G=g|X,S)
parm2 <- ret$CML$parms
if (fixFlag) parm2 <- fixGetEta(parm2)
pmat <- Pdg.xs(Pdg, parm2[alpha.row], parm2[beta.row], parm2[xi.row])$Pdg.matrix
# Get the scores for the 3 parts of eta (parm2)
temp1 <- D - Pd1.xs(pmat)
temp1 <- matrixMultVec(X.main, temp1)
temp2 <- Dsnp - Edg.xs(pmat)
if (gmodel3) {
temp2 <- cbind(matrixMultVec(X.int, temp2[, 1]),
matrixMultVec(X.int, temp2[, 2]))
} else {
temp2 <- matrixMultVec(X.int, temp2)
}
# If the snp was not in the model as a main effect, then remove the
# first column of temp2.
if (zeroSNP) {
if (ncol(temp2) == 1) {
temp2 <- NULL
} else {
temp2 <- temp2[, -1]
}
}
temp3 <- snp - Eg.xs(pmat)
temp3 <- matrixMultVec(X.strata, temp3)
# Get the scores for CML
if (!fit.null) {
score2 <- cbind(temp1, temp2, temp3)
} else {
score2 <- cbind(temp1, temp3)
}
# Free memory
rm(temp1, temp2, temp3, Z, pmat, parm2)
temp <- gc(verbose=FALSE)
# Compute the covariance of beta.UML and eta
temp <- 0
dim1 <- c(ncol(score1), 1)
dim2 <- c(1, ncol(score2))
for (i in 1:n) {
temp1 <- score1[i, ]
temp2 <- score2[i, ]
dim(temp1) <- dim1
dim(temp2) <- dim2
temp <- temp + (temp1 %*% temp2)
}
cov12 <- cov1 %*% temp %*% ret$CML$cov
# We only need the cov(beta.UML, beta.CML) submatrix
cov12 <- cov12[ids, ids]
rm(score1, score2)
temp <- gc(verbose=FALSE)
# Initialize the covariance matrix
nb <- length(parms)
nb2 <- 2*nb
nbp1 <- nb + 1
cov <- matrix(data=NA, nrow=nb2, ncol=nb2)
cov[1:nb, 1:nb] <- cov1
cov[nbp1:nb2, nbp1:nb2] <- ret$CML$cov[ids, ids]
cov[1:nb, nbp1:nb2] <- cov12
cov[nbp1:nb2, 1:nb] <- t(cov12)
# Return this matrix
UML.CML.cov <- cov12
rownames(UML.CML.cov) <- paste("UML.", vnames, sep="")
colnames(UML.CML.cov) <- paste("CML.", vnames, sep="")
# Obtain the final covariance matrix
cov <- cmat %*% cov %*% t(cmat)
colnames(cov) <- vnames
rownames(cov) <- vnames
list(parms=parms, cov=cov, UML.CML.cov=UML.CML.cov)
} # END: getEB
# Function to call before returning the return list
setReturn <- function(ret) {
class(ret) <- "snp.logistic"
if (!is.null(ret$UML)) {
class(ret$UML) <- "UML"
} else {
return(ret)
}
if (is.null(ret$CML)) return(ret)
if (fixFlag) {
xi.row <- fix_xi.row
if (!xi.row[1]) {
class(ret$CML) <- "CML"
return(ret)
}
}
# Transfrom the strata parms ???
xi <- ret$CML$parms[xi.row]
#exi <- exp(xi)
#ret$CML$strata.parms <- exi/(1 + exi)
ret$CML$strata.parms <- xi
# Remove the strata parms from parms
ret$CML$parms <- ret$CML$parms[-xi.row]
# Get the cov matrix for the strata
xi <- ret$CML$cov[xi.row, xi.row]
dim(xi) <- c(nStrata, nStrata)
# Use the delta method. The derivative is a diagonal matrix
#dexi <- ret$CML$strata.parms/(1 + exi)
#ndexi <- length(dexi)
#if (ndexi > 1) dexi <- diag(dexi)
# Get the covariance matrix
#temp <- dexi %*% xi %*% dexi
#dim(temp) <- c(ndexi, ndexi)
#ret$CML$strata.cov <- temp
ret$CML$strata.cov <- xi
# Set the names
temp <- names(ret$CML$strata.parms)
colnames(ret$CML$strata.cov) <- temp
rownames(ret$CML$strata.cov) <- temp
# Save full cov
ret$CML$cov.full <- ret$CML$cov
# Remove strata parms from cov
temp <- ret$CML$cov[-xi.row, -xi.row]
if (length(temp) == 1) {
dim(temp) <- c(1, 1)
rownames(temp) <- colnames(temp) <- "Intercept"
}
ret$CML$cov <- temp
class(ret$CML) <- "CML"
if (!is.null(ret$EB)) class(ret$EB) <- "EB"
ret
} # END: setReturn
# Function to call glm for a factor snp (for now)
fitGLM <- function() {
snp <- factor(snp)
temp <- callGLM(D, X.main=X.main, X.int=X.int, int.vec=snp,
family=op$family, prefix=op$snpName)
if (temp$converged) {
uml <- list(parms=temp$coefficients)
uml$cov <- summary(temp)$cov.scaled
uml$loglike <- getLoglike.glm(temp)
}
list(UML=uml)
} # END: fitGLM
# Function to get a new eta from fixEta0
fixGetEta <- function(eta) {
ret <- fixEta0
ret[fixMap] <- eta
ret
} # END: fixGetEta
# Wrapper for the C function
CML_EB.R <- function() {
if (op$use.C.code == 0) return(NULL)
if (!is.loaded("CML_EB")) {
warning("CML_EB function is not loaded")
return(NULL)
}
# If initial estimates were passed in, then use them
op_eta0 <- op[["init.parms", exact=TRUE]]
if (!is.null(op_eta0)) {
if (length(op_eta0) != nparms) stop("ERROR: init.parms has the wrong length")
eta0[] <- op_eta0
}
error <- 1
# Make the matrices vectors by row
if (!nx) X.main <- 0
if (!nv) X.int <- 0
X.main <- t(X.main)
dim(X.main) <- NULL
X.int <- t(X.int)
dim(X.int) <- NULL
X.strata <- t(X.strata)
dim(X.strata) <- NULL
# Define the return vector
cml.parms <- double(nparms)
cml.cov <- double(nparms*nparms)
cml.ll <- double(1)
error <- as.integer(1)
nbp1 <- nbeta + 1
eb.parms <- double(nbp1)
eb.cov <- double(nbp1*nbp1)
uml.cov <- ret$UML$cov
uml.parms <- ret$UML$parms
if (zeroSNP) {
# Remove snp and update other variables
temp <- match(op$snpName, names(eta0))
if (!is.na(temp)) {
eta0 <- eta0[-temp]
nparms <- length(eta0)
nbeta <- nbeta - 1
}
}
# Vector for UML-CML cov
uml.cml.cov <- double((nbeta+1)*nparms)
##########################################################################
############## Include PACKAGE="CGEN" when building a package ############
##########################################################################
# Call the C function
temp <- .C("CML_EB", as.double(eta0), as.integer(nparms), as.integer(nbeta),
as.integer(D), as.double(snp), as.integer(n), as.double(X.main), as.integer(nx),
as.double(X.int), as.integer(nv), as.double(X.strata), as.integer(nStrata),
as.integer(genetic.model), as.integer(geno.binary),
as.integer(op$maxiter), as.double(op$reltol), as.integer(op$debug),
as.double(uml.cov), as.double(fitted.values),
as.integer(zeroSNP), as.integer(op$num.deriv), as.integer(imputed), as.double(ProbG1),
as.double(uml.parms), error=error, cml.parms=cml.parms, cml.cov=cml.cov, cml.ll=cml.ll,
eb.parms=eb.parms, eb.cov=eb.cov, uml.cml.cov=uml.cml.cov, PACKAGE="CGEN")
error <- temp$error
if (error) return(list())
# Get the covariance matrix
cml.cov <- matrix(temp$cml.cov, nrow=nparms, byrow=TRUE)
if (any(!is.finite(diag(cml.cov)))) return(list())
cml.parms <- temp$cml.parms
cnames <- names(eta0)
names(cml.parms) <- cnames
colnames(cml.cov) <- cnames
rownames(cml.cov) <- cnames
cnames <- cnames[1:(1+nbeta)]
eb.parms <- temp$eb.parms
names(eb.parms) <- cnames
eb.cov <- matrix(temp$eb.cov, nrow=1+nbeta, byrow=TRUE)
colnames(eb.cov) <- cnames
rownames(eb.cov) <- cnames
# UML-CML matrix
uml.cml.cov <- matrix(temp$uml.cml.cov, byrow=TRUE, ncol=nparms)
rownames(uml.cml.cov) <- paste("UML.", cnames, sep="")
colnames(uml.cml.cov) <- paste("CML.", names(eta0), sep="")
# Return list
list(CML=list(parms=cml.parms, cov=cml.cov, loglike=temp$cml.ll),
EB=list(parms=eb.parms, cov=eb.cov, UML.CML.cov=uml.cml.cov))
} # END: CML_EB.R
# Function to check the initial estimates of CML
check_init <- function() {
new <- eta0
# If initial estimates were passed in, then use them
op_eta0 <- op[["init.parms", exact=TRUE]]
if (!is.null(op_eta0)) {
if (length(op_eta0) != nparms) stop("ERROR: init.parms has the wrong length")
new[] <- op_eta0
}
maxll <- loc.getLoglike(eta0)
maxit <- 100
h <- 0.1
steps <- h*abs(new)
temp <- steps <= 1e-8
steps[temp] <- 0.01
for (i in 1:nparms) {
test <- new
point <- new[i]
step <- steps[i]
flag <- 0
# For 2 directions
for (k in 1:2) {
for (j in 1:maxit) {
point0 <- point + step
test[i] <- point0
ll <- loc.getLoglike(test)
if (ll > maxll) {
maxll <- ll
point <- point0
flag <- 1
} else {
break
}
}
if (flag) {
new[i] <- point
break
} else {
# Try other direction
point <- new[i]
step <- -step
}
}
}
new
} # END: check_init
# Check the options list
op <- default.list(op, c("reltol", "maxiter", "optimizer",
"snpName", "debug", "genetic.model", "errorCheck", "geno.binary",
"use.C.code", "only.UML", "fit.null", "num.deriv", "imputed",
"s.1catVar"),
list(1e-8, 100, "BFGS", "SNP_", 0, 0, 1, 0, 1, 0, 0, 0, 0, 0))
snp.nc <- ncol(snp)
if (is.null(snp.nc)) snp.nc <- 0
if (snp.nc == 3) op$imputed <- 1
op$imputed <- 0 # Changed Mar 11, 2015
fixFlag <- 0
fit.null <- op$fit.null
imputed <- op$imputed
if (imputed) {
if (snp.nc == 3) {
ProbG1 <- snp[, 2]
# Create snp vector
gmodel <- op$genetic.model
if (gmodel == 0) {
snp <- snp[, 2] + 2*snp[, 3]
} else if (gmodel == 1) {
snp <- snp[, 2] + snp[, 3]
} else if (gmodel == 2) {
snp <- snp[, 3]
} else {
stop("Incorrect genetic.model with imputed data")
}
}
op$genetic.model <- 0
if (is.null(ProbG1)) stop("ERROR: ProbG1 is NULL")
} else {
ProbG1 <- 0
}
zeroSNP <- FALSE
if (fit.null) {
X.int <- NULL
op$fixParms <- list(parms=op$snpName, values=0)
zeroSNP <- TRUE
fixFlag <- 1
} else {
fixFlag <- !is.null(op[["fixParms", exact=TRUE]])
if (fixFlag) {
if (op$snpName %in% op$fixParms$parms) zeroSNP <- TRUE
}
}
if (zeroSNP) op$genetic.model <- 0
genetic.model <- op$genetic.model
geno.binary <- 0
if (!(genetic.model %in% c(0, 1, 2, 3))) stop("genetic.model must be 0-3")
# Get the number of genotypes
snp <- unfactor(snp, fun=as.numeric)
usnp <- sort(unique(snp))
n <- length(usnp)
# If the input SNP is binary 0-1, set genetic.model to 0
if (!imputed) {
if (!all(usnp %in% c(0, 1, 2))) stop("snp is not coded correctly")
}
if (n == 1) {
stop("snp only has 1 value")
} else if (n == 2) {
if (genetic.model) warning("genetic.model is set to 0")
genetic.model <- 0
if (all(usnp %in% 0:1)) geno.binary <- 1
}
if (genetic.model %in% 1:2) geno.binary <- 1
# Check D
if (!all(unique(D) %in% c(0, 1))) stop("D is not coded correctly")
if (!is.null(X.strata)) {
if (!is.matrix(X.strata)) stop("X.strata is not a matrix")
}
gmodel3 <- (genetic.model == 3)
n <- length(D)
if (op$optimizer != "BFGS") op$use.C.code <- 0
#if (imputed) dim(ProbG1) <- c(n, 1)
# See if X and V were given
if (is.null(X.main)) {
nx <- 0
} else {
nx <- ncol(X.main)
}
if (is.null(X.int)) {
nv <- 0
} else {
nv <- ncol(X.int)
}
# Get the SNP vector for the specific genetic model
fsnp <- snp
if (genetic.model == 1) {
# Dominant
temp <- snp == 2
fsnp[temp] <- 1
} else if (genetic.model == 2) {
temp <- snp == 1
fsnp[temp] <- 0
temp <- snp == 2
fsnp[temp] <- 1
} else if (genetic.model == 3) {
fsnp <- cbind(as.integer(snp == 1), as.integer(snp == 2))
}
# Check the strata
if (is.null(X.strata)) X.strata <- matrix(data=1, nrow=n, ncol=1)
nStrata <- ncol(X.strata)
# Get the initial estimates
temp <- getInit()
if (is.null(temp)) return(NULL)
# Save the initial estimates and initialize the return list
ret <- list()
ret$UML <- list(parms=temp$UML.parms, cov=temp$cov, loglike=temp$loglike)
if (op$only.UML) setReturn(ret)
nbeta <- length(temp$beta)
eta0 <- temp$eta
alpha.row <- 1
beta.row <- 2:(nbeta+1)
xi.row <- (max(beta.row)+1):length(eta0)
cov <- NULL
nparms <- length(eta0)
# Save the fitted values for empirical bayes
fitted.values <- temp$fitted.values
# If there were NAs in UML, modify the design matrices
if (temp$naFlag) {
if (nx) {
pos <- temp$naXPos
len <- length(pos)
if (len) {
if (len == nx) {
nx <- 0
X.main <- NULL
} else {
X.main <- removeOrKeepCols(X.main, pos, which=-1)
nx <- ncol(X.main)
}
}
}
if (nv) {
pos <- temp$naVPos
len <- length(pos)
if (len) {
if (len == nv) {
nv <- 0
X.int <- NULL
} else {
X.int <- removeOrKeepCols(X.int, pos, which=-1)
nv <- ncol(X.int)
}
}
}
}
# Determine if parameters are to be fixed
if (fixFlag) {
#op$use.C.code <- 0
temp <- op$fixParms
# Define the map for fixed parms
temp2 <- match(temp$parms, names(eta0))
temp2 <- temp2[!is.na(temp2)]
if (!length(temp2)) stop("ERROR with fixParms$parms")
fixEta0 <- eta0
fixEta0[temp2] <- temp$values
fixMap <- (1:nparms)[-temp2]
# Change eta0
eta0 <- eta0[fixMap]
# Get the updated xi.row
temp <- sum(xi.row %in% temp2)
nxi <- length(xi.row) - temp
if (nxi) {
temp <- length(eta0)
fix_xi.row <- (temp-nxi+1):temp
} else {
fix_xi.row <- 0
}
}
if (genetic.model) op$num.deriv <- 1
# Call the C code
clist <- try(CML_EB.R(), silent=TRUE)
if (checkTryError(clist, conv=0)) return(setReturn(ret))
if (!is.null(clist)) {
if (!length(clist)) return(setReturn(ret))
ret$CML <- clist$CML
ret$EB <- clist$EB
return(setReturn(ret))
}
# Initialize
nlevels <- 6
d0g0.col <- 1
d0g1.col <- 2
d0g2.col <- 3
d1g0.col <- 4
d1g1.col <- 5
d1g2.col <- 6
d0.col <- c(d0g0.col, d0g1.col, d0g2.col)
d1.col <- c(d1g0.col, d1g1.col, d1g2.col)
g0.col <- c(d0g0.col, d1g0.col)
g1.col <- c(d0g1.col, d1g1.col)
g2.col <- c(d0g2.col, d1g2.col)
log2 <- log(2)
# Define the Z matrices for efficiency
Z0 <- getZ(0)
Z1 <- getZ(1)
Z2 <- getZ(2)
if (imputed) Z.imp <- getZ.imp(fsnp)
# Initialize the matrix to hold all probabilities P(D=d, G=g| X, S)
Pdg <- matrix(data=0, nrow=n, ncol=nlevels)
# Compute D*snp
if (gmodel3) {
Dsnp <- matrixMultVec(fsnp, D)
} else {
Dsnp <- D*fsnp
}
# Add intercept columns to X and V
X.main <- addIntercept(X.main, nrow=n)
X.int <- addIntercept(X.int, nrow=n)
# Define a logical matrix for the calculation of the log-likelihood
loglike.mat <- getLoglike.mat()
# Update initial estimates if needed
eta0 <- check_init()
temp <- try(callOptim(), silent=TRUE)
if (checkTryError(temp, conv=0)) return(setReturn(ret))
if (!temp$converged) return(setReturn(ret))
temp <- list(parms=temp$parms, cov=temp$cov, loglike=temp$loglike)
ret$CML <- temp
# Empirical Bayes
ret$EB <- getEB()
ret <- setReturn(ret)
ret
} # END: snp.main
# Function to return a vector of variable names
logistic.vnames <- function(X, V, nStrata, snpName="SNP_",
out.est=NULL, genetic.model=0,
fit.null=0, zeroSNP=0) {
# Initialize the return list
ret <- list()
# Check out.est
if (!is.null(out.est)) {
temp <- out.est$parms
if (is.character(temp)) {
ret$est.p <- temp
out.est <- 0
} else {
# Numeric value
out.est <- temp
}
} else {
out.est <- 0
}
save1 <- out.est
ret$int <- "Intercept"
if (!is.null(X)) ret$X <- getVarNames(X, prefix="X")
ret$snp <- getVarNames.snp(prefix=snpName, genetic.model=genetic.model)
if (!is.null(V)) {
temp <- getVarNames.int(V, prefix=snpName, genetic.model=genetic.model, sep=":")
} else {
temp <- NULL
if (out.est == 2) out.est <- 1
}
ret$V <- temp
ret$strata <- paste("Allele_freq.", 1:nStrata, sep="")
ret$all <- c(ret$int, ret$X, ret$snp, ret$V, ret$strata)
if (out.est) {
if (out.est == 3) {
ret$est.p <- c(ret$int, ret$X, ret$snp, ret$V)
} else {
ret$est.p <- ret$snp
if (out.est == 2) {
ret$est.p <- c(ret$est.p, ret$V)
}
}
if (genetic.model == 3) {
# Create another vector of parm names for the case when there is
# less than 3 genotypes
snp1 <- paste(snpName, 1, sep="")
if (!is.null(V)) {
tempv <- getVarNames.int(V, prefix=snpName, genetic.model=0, sep=":")
tempv1 <- getVarNames.int(V, prefix=snp1, genetic.model=0, sep=":")
} else {
tempv <- NULL
tempv1 <- NULL
}
if (save1 == 3) {
temp <- c(ret$int, ret$X, snpName, tempv)
temp1 <- c(ret$int, ret$X, snp1, tempv1)
} else if (save1 == 2) {
temp <- c(snpName, tempv)
temp1 <- c(snp1, tempv1)
} else if (save1 == 1) {
temp <- snpName
temp1 <- snp1
}
ret$est.p <- list(ret$est.p, temp)
ret$est.parms <- list(ret$est.p[[1]], temp1)
} else {
ret$est.parms <- ret$est.p
}
} else {
ret$est.parms <- ret$est.p
}
if (fit.null) {
ret$UML.names <- c(ret$int, ret$X)
ret$all.names <- c(ret$int, ret$X, ret$snp, ret$strata)
} else {
UML.names <- c(ret$int, ret$X)
if (!zeroSNP) UML.names <- c(UML.names, ret$snp)
UML.names <- c(UML.names, ret$V)
ret$UML.names <- UML.names
ret$all.names <- ret$all
}
ret
} # END: logistic.vnames
# Function to create a design matrix
logistic.dsgnMat <- function(data, vars, facVars, removeInt=1,
norm.names=1, remove.vars=NULL) {
# data Data frame
# vars Character vector of variable names or a formula
# facVars Character vector of factor names
# remove.vars Character vector of variable names to remove.
# They are removed after the variable names are normalized.
# The default is NULL.
rm.vars <- function(mat, rmvars) {
mnames <- colnames(mat)
temp <- !(mnames %in% rmvars)
mnames <- mnames[temp]
len <- length(mnames)
if (len == 0) {
mat <- NULL
} else if ((len == 1) && (mnames == "Intercept")) {
mat <- NULL
} else {
mat <- removeOrKeepCols(mat, mnames, which=1)
}
mat
} # END: rm.vars
if (is.null(vars)) return(list(designMatrix=NULL, newVars=NULL))
# Determine if vars is a formula
if ("formula" %in% class(vars)) {
# Get the design matrix
design <- model.matrix(vars, data=data)
# Remove the intercept, if needed
newVars <- colnames(design)
if (removeInt) {
if (newVars[1] == "(Intercept)") {
design <- removeOrKeepCols(design, 1, which=-1)
newVars <- newVars[-1]
}
}
if (norm.names) colnames(design) <- normVarNames(colnames(design))
if (!is.null(remove.vars)) design <- rm.vars(design, remove.vars)
return(list(designMatrix=design, newVars=newVars))
}
design <- removeOrKeepCols(data, vars, which=1)
newVars <- NULL
if (!is.null(facVars)) {
temp <- vars %in% facVars
if (any(temp)) {
temp <- vars[temp]
temp <- createDummy(design, vars=temp)
design <- temp$data
newVars <- temp$newVars
}
}
design <- as.matrix(design)
if (norm.names) colnames(design) <- normVarNames(colnames(design))
# Check for constant variables
design <- checkForConstantVar(design, msg=1)$data
if (!removeInt) {
# Add intercept
cnames <- colnames(design)
design <- cbind(1, design)
colnames(design) <- c("Intercept", cnames)
}
if (!is.null(remove.vars)) design <- rm.vars(design, remove.vars)
# Make sure matrix is numeric
d <- dim(design)
cnames <- colnames(design)
design <- as.numeric(design)
dim(design) <- d
colnames(design) <- cnames
list(designMatrix=design, newVars=newVars)
} # END: logistic.dsgnMat
# Function to apply the zero.vars option to the design matrices
apply_zero.vars <- function(zero.vars, mat.list, snp.var, facVars, data) {
mat.names <- names(mat.list)
fixParms <- NULL
if (is.character(zero.vars)) {
for (i in 1:length(mat.names)) {
mat <- mat.list[[i]]
if (is.null(mat)) next
cnames <- colnames(mat)
temp <- cnames %in% zero.vars
if (any(temp)) {
cnames <- cnames[!temp]
if (!length(cnames)) {
mat <- NULL
} else {
mat <- removeOrKeepCols(mat, cnames, which=1)
}
mat.list[[mat.names[i]]] <- mat
}
}
if (snp.var %in% zero.vars) {
fixParms <- list(parms=snp.var, values=0)
}
} else {
# zero.vars is a list
znames <- names(zero.vars)
temp <- zero.vars[["snp.var", exact=TRUE]]
if (!is.null(temp)) {
if (temp == snp.var) {
fixParms <- list(parms=snp.var, values=0)
}
}
for (z in znames) {
vars <- zero.vars[[z]]
if (is.null(vars)) next
mat <- mat.list[[z, exact=TRUE]]
if (is.null(mat)) next
cnames <- colnames(mat)
pnames <- colnames(logistic.dsgnMat(data, vars, facVars)$designMatrix)
temp <- cnames %in% pnames
if (any(temp)) {
cnames <- cnames[!temp]
if (!length(cnames)) {
mat <- NULL
} else {
mat <- removeOrKeepCols(mat, cnames, which=1)
}
mat.list[[z]] <- mat
}
}
}
list(mat.list=mat.list, fixParms=fixParms)
} # END: apply_zero.vars
# Function to output the needed UML and CML estimates
UML_CML_GxE_parms <- function(main.variables, interaction.variables, out.file, snps, data, response.var,
main.vars=NULL, int.vars=NULL, strata.var=NULL, ProbG1.var=NULL, op=NULL) {
# Input arguments:
# out.file Output file to save the necessary UML and CML estimates, or NULL
# snps Character vector of snps to be analyzed
# data Data frame containing the snps to be analyzed, outcome,
# covariates and other variables
# response.var Same definition as is snp.logistic
# main.vars Same definition as is snp.logistic
# int.var Character vector of interaction variable
# strata.var Same definition as is snp.logistic
# op Same definition as is snp.logistic
# Output variable names for the UML main effects:
# UML.G.BETA UML.E.BETA UML.GE.BETA
# Output variable names for the upper triangle of the UML covariance matrix:
# UML.G.G.COV, UML.G.E.COV UML.G.GE.COV
# UML.E.E.COV UML.E.GE.COV
# UML.GE.GE.COV
# Output variable names for the UML-CML covariance matrix:
# UML.G.CML.G.COV, UML.G.CML.E.COV UML.G.CML.GE.COV
# UML.E.CML.G.COV, UML.E.CML.E.COV UML.E.CML.GE.COV
# UML.GE.CML.G.COV, UML.GE.CML.E.COV UML.GE.CML.GE.COV
n.main.vars <- length(main.variables)
n.int.vars <- length(interaction.variables)
UML.ME.out <- outNames.me("UML", n.main.vars, n.int.vars)
CML.ME.out <- outNames.me("CML", n.main.vars, n.int.vars)
UML.COV.out <- outNames.cov("UML", n.main.vars, n.int.vars)
CML.COV.out <- outNames.cov("CML", n.main.vars, n.int.vars)
UML.CML.COV.out <- outNames.cov2(n.main.vars, n.int.vars)
fileFlag <- !is.null(out.file)
print <- op$print
# Open the output file
if (fileFlag) fid <- file(out.file, "w")
# Write out the column names
outvars <- c("SNP", UML.ME.out, UML.COV.out, CML.ME.out, CML.COV.out, UML.CML.COV.out)
if (fileFlag) writeOut(fid, outvars)
# Vector to store output values
outvec <- character(length(outvars))
names(outvec) <- outvars
outvec0 <- outvec
# Loop over each snp
for (snp in snps) {
outvec0[] <- ""
outvec <- outvec0
outvec["SNP"] <- snp
fit <- try(snp.logistic(data, response.var, snp, main.vars=main.vars,
int.vars=int.vars, strata.var=strata.var, op=op), silent=TRUE)
if ("try-error" %in% class(fit)) {
if (print) print(fit)
if (fileFlag) writeOut(fid, outvec)
next
}
if (print) print(summary(fit))
# Extract UML and CML estimates
temp <- try(extractEst(fit, "UML", outvec, snp, main.variables, interaction.variables,
UML.ME.out, UML.COV.out), silent=TRUE)
if (!("try-error" %in% class(temp))) outvec <- temp
temp <- try(extractEst(fit, "CML", outvec, snp, main.variables, interaction.variables,
CML.ME.out, CML.COV.out), silent=TRUE)
if (!("try-error" %in% class(temp))) outvec <- temp
temp <- try(extract_UML.CML(fit, outvec, snp, main.variables, interaction.variables,
UML.CML.COV.out), silent=TRUE)
if (!("try-error" %in% class(temp))) outvec <- temp
if (fileFlag) writeOut(fid, outvec)
}
if (fileFlag) close(fid)
outvec
} # END: UML_CML_GxE_parms
# Function to write output to an open (tab-delimited) file
writeOut <- function(FID, vec) {
# FID File id
# vec Vector to output
str <- paste(vec, collapse="\t", sep="")
write(str, file=FID, ncolumns=1)
NULL
} # END: writeOut
# Function to get the parameter names
outNames.me <- function(which, n.main.vars, n.int.vars) {
ret <- "G"
if (n.main.vars < 2) {
ret <- c(ret, "E")
} else {
jj <- 1:n.main.vars
ret <- c(ret, paste("E", jj, sep=""))
}
if (n.int.vars < 2) {
ret <- c(ret, "GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("GE", jj, sep=""))
}
ret <- paste(which, ".", ret, ".BETA", sep="")
ret
} # END: outNames.me
# Function to get output names
outNames.cov <- function(which, n.main.vars, n.int.vars) {
# G row
ret <- "G.G"
if (n.main.vars < 2) {
ret <- c(ret, "G.E")
} else {
jj <- 1:n.main.vars
ret <- c(ret, paste("G.E", jj, sep=""))
}
if (n.int.vars < 2) {
ret <- c(ret, "G.GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("G.GE", jj, sep=""))
}
# E rows
if (n.main.vars < 2) {
ret <- c(ret, "E.E")
if (n.int.vars < 2) {
ret <- c(ret, "E.GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("E.GE", jj, sep=""))
}
} else {
jj0 <- 1:n.main.vars
jj2 <- 1:n.int.vars
for (i in 1:n.main.vars) {
jj <- i:n.main.vars
ret <- c(ret, paste("E", i, ".E", jj, sep=""))
if (n.int.vars < 2) {
ret <- c(ret, paste("E", i, ".GE", sep=""))
} else {
ret <- c(ret, paste("E", i, ".GE", jj2, sep=""))
}
}
}
# GE rows
if (n.int.vars < 2) {
ret <- c(ret, "GE.GE")
} else {
jj0 <- 1:n.int.vars
for (i in 1:n.int.vars) {
jj <- i:n.int.vars
ret <- c(ret, paste("GE", i, ".GE", jj, sep=""))
}
}
ret <- paste(which, ".", ret, ".COV", sep="")
ret
} # END: outNames.cov
outNames.cov2 <- function(n.main.vars, n.int.vars) {
if ((n.main.vars < 2) && (n.int.vars < 2)) {
ret <- c("UML.G.CML.G.COV", "UML.G.CML.E.COV", "UML.G.CML.GE.COV",
"UML.E.CML.G.COV", "UML.E.CML.E.COV", "UML.E.CML.GE.COV",
"UML.GE.CML.G.COV", "UML.GE.CML.E.COV", "UML.GE.CML.GE.COV")
return(ret)
}
# G row
ret <- "UML.G.CML.G"
if (n.main.vars < 2) {
ret <- c(ret, "UML.G.CML.E")
} else {
jj <- 1:n.main.vars
ret <- c(ret, paste("UML.G.CML.E", jj, sep=""))
}
if (n.int.vars < 2) {
ret <- c(ret, "UML.G.CML.GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("UML.G.CML.GE", jj, sep=""))
}
# E rows
if (n.main.vars < 2) {
ret <- c(ret, "UML.E.CML.G")
ret <- c(ret, "UML.E.CML.E")
if (n.int.vars < 2) {
ret <- c(ret, "UML.E.CML.GE")
} else {
ret <- c(ret, paste("UML.E.CML.GE", 1:n.int.vars, sep=""))
}
} else {
jj <- 1:n.main.vars
jj2 <- 1:n.int.vars
for (i in 1:n.main.vars) {
ret <- c(ret, paste("UML.E", i, ".CML.G", sep=""))
ret <- c(ret, paste("UML.E", i, ".CML.E", jj, sep=""))
if (n.int.vars < 2) {
ret <- c(ret, paste("UML.E", i, ".CML.GE", sep=""))
} else {
ret <- c(ret, paste("UML.E", i, ".CML.GE", jj2, sep=""))
}
}
}
# GE rows
if (n.int.vars < 2) {
ret <- c(ret, "UML.GE.CML.G")
if (n.main.vars < 2) {
ret <- c(ret, "UML.GE.CML.E")
} else {
ret <- c(ret, paste("UML.GE.CML.E", 1:n.main.vars, sep=""))
}
ret <- c(ret, "UML.GE.CML.GE")
} else {
jj0 <- 1:n.main.vars
jj2 <- 1:n.int.vars
for (i in 1:n.int.vars) {
ret <- c(ret, paste("UML.GE", i, ".CML.G", sep=""))
if (n.main.vars < 2) {
ret <- c(ret, paste("UML.GE", i, ".CML.E", sep=""))
} else {
ret <- c(ret, paste("UML.GE", i, ".CML.E", jj0, sep=""))
}
ret <- c(ret, paste("UML.GE", i, ".CML.GE", jj2, sep=""))
}
}
ret <- paste(ret, ".COV", sep="")
ret
} # END: outNames.cov2
# Function to get the UML-CML cov estimates
extract_UML.CML <- function(fit, outvec, G.name, E.name, GE.name, ids) {
x <- fit[["EB", exact=TRUE]]
if (is.null(x)) return(outvec)
cov <- x[["UML.CML.cov", exact=TRUE]]
if (is.null(cov)) return(outvec)
# Get the names of the interaction parms
GE.name <- paste(G.name, ":", GE.name, sep="")
G.row <- paste("UML.", G.name, sep="")
G.col <- paste("CML.", G.name, sep="")
E.row <- paste("UML.", E.name, sep="")
E.col <- paste("CML.", E.name, sep="")
GE.row <- paste("UML.", GE.name, sep="")
GE.col <- paste("CML.", GE.name, sep="")
# Vector of expected parameter names
row.names <- paste("UML.", c(G.name, E.name, GE.name), sep="")
col.names <- paste("CML.", c(G.name, E.name, GE.name), sep="")
n.names <- length(col.names)
n.int.vars <- length(E.name)
cnames <- colnames(cov)
temp <- col.names %in% cnames
if (!all(temp)) {
cov2 <- matrix(data=NA, nrow=n.names, ncol=n.names)
colnames(cov2) <- col.names
rownames(cov2) <- row.names
col2.names <- col.names[temp]
row2.names <- row.names[temp]
cov2[row2.names, col2.names] <- cov[row2.names, col2.names]
} else {
cov2 <- cov
}
if (n.int.vars < 2) {
outvec[ids] <- c(cov2[G.row, G.col], cov2[G.row, E.col], cov2[G.row, GE.col],
cov2[E.row, G.col], cov2[E.row, E.col], cov2[E.row, GE.col],
cov2[GE.row, G.col], cov2[GE.row, E.col], cov2[GE.row, GE.col])
} else {
vec <- cov2[G.row, G.col]
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.row, E.col[i]])
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.row, GE.col[i]])
for (i in 1:n.int.vars) {
row <- E.row[i]
vec <- c(vec, cov2[row, G.col])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, E.col[j]])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, GE.col[j]])
}
for (i in 1:n.int.vars) {
row <- GE.row[i]
vec <- c(vec, cov2[row, G.col])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, E.col[j]])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, GE.col[j]])
}
outvec[ids] <- vec
}
outvec
} # END: extract_UML.CML
# Function to extract estimates. The updated output vector will be returned
extractEst <- function(fit, which, outvec, G.name, E.name, GE.name, ids.main, ids.cov) {
# fit Return object from snp.logistic
# which "UML" or "CML"
# outvec Output vector
# G.name Name of the SNP
# E.name Name of the environmental variable
if (which == "EB") return(extract_UML.CML(fit, outvec, G.name, E.name))
x <- fit[[which, exact=TRUE]]
if (is.null(x)) return(outvec)
n.int.vars <- length(E.name)
# Get the names of the interaction parms
GE.name <- paste(G.name, ":", GE.name, sep="")
# Vector of expected parameter names
vec.names <- c(G.name, E.name, GE.name)
n.vec.names <- length(vec.names)
parms <- x[["parms", exact=TRUE]]
if (!is.null(parms)) {
cnames <- names(parms)
temp <- vec.names %in% cnames
if (!all(temp)) {
parms2 <- rep(NA, n.vec.names)
names(parms2) <- vec.names
vec2.names <- vec.names[temp]
parms2[vec2.names] <- parms[vec2.names]
} else {
parms2 <- parms
vec2.names <- vec.names
}
outvec[ids.main] <- c(parms[vec2.names])
}
cov <- x[["cov", exact=TRUE]]
if (!is.null(cov)) {
cnames <- colnames(cov)
temp <- vec.names %in% cnames
if (!all(temp)) {
cov2 <- matrix(data=NA, nrow=n.vec.names, ncol=n.vec.names)
colnames(cov2) <- vec.names
rownames(cov2) <- vec.names
vec2.names <- vec.names[temp]
cov2[vec2.names, vec2.names] <- cov[vec2.names, vec2.names]
} else {
cov2 <- cov
}
if (n.int.vars < 2) {
outvec[ids.cov] <- c(cov2[G.name, G.name], cov2[G.name, E.name], cov2[G.name, GE.name],
cov2[E.name, E.name], cov2[E.name, GE.name], cov2[GE.name, GE.name])
} else {
vec <- cov2[G.name, G.name]
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.name, E.name[i]])
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.name, GE.name[i]])
for (i in 1:n.int.vars) {
row <- E.name[i]
for (j in i:n.int.vars) vec <- c(vec, cov2[row, E.name[j]])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, GE.name[j]])
}
for (i in 1:n.int.vars) {
row <- GE.name[i]
for (j in i:n.int.vars) vec <- c(vec, cov2[row, GE.name[j]])
}
outvec[ids.cov] <- vec
}
}
outvec
} # END: extractEst
# Function for running a scan
scan.UML_CML <- function(snp.list, pheno.list, op=NULL) {
# Check the input lists
snp.list <- check.snp.list(snp.list)
pheno.list <- default.list(pheno.list,
c("response.var"), list("ERROR"), error=c(1))
response.var <- pheno.list$response.var
main.vars <- pheno.list[["main.vars", exact=TRUE]]
int.vars <- pheno.list[["int.vars", exact=TRUE]]
strata.var <- pheno.list[["strata.vars", exact=TRUE]]
cc.var <- pheno.list[["cc.var", exact=TRUE]]
# Determine if fomulas were passed in
temp <- main.vars
if (!is.null(temp)) {
main.form <- ("formula" %in% class(temp))
} else {
main.form <- FALSE
}
temp <- int.vars
if (!is.null(temp)) {
int.form <- ("formula" %in% class(temp))
} else {
int.form <- FALSE
}
temp <- strata.var
if (!is.null(temp)) {
s.form <- ("formula" %in% class(temp))
} else {
s.form <- FALSE
}
formFlag <- main.form + int.form + s.form
temp.list <- op[["temp.list", exact=TRUE]]
temp.list <- check.temp.list(temp.list)
temp <- c("BFGS", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "IRLS")
op <- default.list(op,
c("id", "print", "optimizer", "snpName","genetic.model",
"allele.cc", "out.file", "save.varnames"),
list(1, 0, "BFGS", "SNP_", 0, 1, "scan.UML_CML.txt", 1),
error=c(0, 0, 0, 0, 0, 0, 0, 0),
checkList=list(NA, 0:1, temp, NA, 0:2, 1:2, NA, 0:1))
print <- op$print
SNP <- op$snpName
# Check cc.var
if (is.null(cc.var)) {
if (op$allele.cc == 2) pheno.list$cc.var <- pheno.list$response.var
}
# The genetic model is taken care of in logistic.SNP
snp.list$genetic.model <- NULL
# All the variables must be given by variable name (not column number)
if ((is.numeric(pheno.list$response.var)) ||
(is.numeric(pheno.list$strata.var)) ||
(is.numeric(pheno.list$factor.vars)) ||
(is.numeric(pheno.list$id.var)) ) {
stop("ERROR: variables must be specified by name, not column number")
}
if ((!main.form) && (is.numeric(pheno.list$X.vars)))
stop("ERROR: variables must be specified by name, not column number")
if ((!int.form) && (is.numeric(pheno.list$V.vars)))
stop("ERROR: variables must be specified by name, not column number")
# Keep only the variables we need
if (!formFlag) {
temp <- c(response.var, strata.var, main.vars, int.vars,
pheno.list$id.var, cc.var)
pheno.list$keep.vars <- unique(temp)
pheno.list$remove.vars <- NULL
} else {
# Do not remove variables
pheno.list$keep.vars <- NULL
}
pheno.list$remove.miss <- 1
pheno.list$make.dummy <- 0
# Get the data vector of snps
tlist <- list(include.row1=0, include.snps=0, return.type=1, MAF=1,
missing=1, snpNames=1, orderByPheno=1, return.pheno=1, ProbG1=1)
temp <- try(getData.1(snp.list, pheno.list, temp.list, op=tlist),
silent=TRUE)
if (class(temp) == "try-error") {
print(temp)
stop("ERROR loading data")
}
snpData <- temp$data
snpNames <- temp$snpNames
delimiter <- getDelimiter(snp.list, output=1)
nsnps <- length(snpData)
maf.vec <- temp[["MAF", exact=TRUE]]
maf.flag <- !is.null(maf.vec)
alleles <- temp[["alleles", exact=TRUE]]
allFlag <- !is.null(alleles)
ProbG1 <- temp[["ProbG1", exact=TRUE]]
ProbG1.flag <- !is.null(ProbG1)
ProbG1.var <- NULL
# Get the phenotype data
phenoData.list <- temp$phenoData.list
phenoData0 <- phenoData.list$data
phenoData0 <- as.data.frame(phenoData0, stringsAsFactors=FALSE)
phenoData0[, SNP] <- NA
if (ProbG1.flag) {
ProbG1.var <- paste(SNP, "ProbG1", sep="")
phenoData0[, ProbG1.var] <- NA
}
rm(temp, tlist, snp.list, phenoData.list)
temp <- gc(verbose=FALSE)
# Get the response variable
response0 <- as.numeric(phenoData0[, response.var])
if (!any(response0 %in% 0:1)) stop("ERROR: response must be 0-1")
# Get the number of observations
nobs <- length(response0)
# Print out number of observations that will be used
temp <- paste("For the analysis, ", nobs, " observations will be used.\n", sep="")
cat(temp)
cat(paste("Number of cases = ", sum(response0 %in% 1), "\n", sep=""))
cat(paste("Number of controls = ", sum(response0 %in% 0), "\n", sep=""))
cat(paste("Output will be written to: ", op$out.file, "\n\n", sep=""))
rm(pheno.list)
temp <- gc(verbose=FALSE)
# Run a base model with simulated SNP
phenoData0[, SNP] <- rbinom(nobs, 1, 0.5)
out.file <- op$out.file
out.base <- paste(out.file, "_info", sep="")
temp <- try(snp.logistic(phenoData0, response.var, SNP, main.vars=main.vars,
int.vars=int.vars, strata.var=NULL, op=op), silent=TRUE)
if ("try-error" %in% class(temp)) {
print(temp)
stop()
}
interaction.variables <- temp$model.info$int.vars
main.variables <- op[["E.parm.names", exact=TRUE]]
if (is.null(main.variables)) main.variables <- interaction.variables
if (op$save.varnames) {
sink(out.base)
cat("Exposure variables:\n")
print(main.variables)
cat("Interaction variables:\n")
print(interaction.variables)
print(summary(temp))
print(table(response0))
for (v in unique(c(main.variables, interaction.variables))) {
cat(paste(v, "\n", sep=""))
tab <- table(phenoData0[, v], response0, exclude=NULL)
if (length(tab) < 31) {
tab <- addmargins(tab)
print(tab)
}
}
sink()
}
if (print) {
cat("Exposure variables:\n")
print(main.variables)
cat("Interaction variables:\n")
print(interaction.variables)
}
rm(response0)
gc()
# Make sure all interaction variables are also main effects
temp <- temp$UML$parms
temp <- names(temp)
if (!(all(main.variables %in% temp))) {
stop("All interaction variables must also be main effects")
}
# Flag for first line of output
firstFlag <- 0
# Open the output file
fid <- file(op$out.file, "w")
# Loop over each SNP
i <- 0
while (1) {
errorFlag <- 0
i <- i + 1
if (i > nsnps) break
snp <- as.numeric(getVecFromStr(snpData[i], delimiter=delimiter))
snp.name <- snpNames[i]
phenoData0[, SNP] <- snp
if (allFlag) {
majMin <- alleles[i]
} else {
majMin <- " "
}
# Get the MAF
if (maf.flag) {
maf <- maf.vec[i]
} else {
maf <- getMAF(snp)
}
# For imputed data
if (ProbG1.flag) {
temp <- as.numeric(getVecFromStr(ProbG1[i], delimiter=delimiter))
phenoData0[, ProbG1.var] <- temp
}
# Fit the model
temp <- try(UML_CML_GxE_parms(main.variables, interaction.variables, NULL, SNP, phenoData0, response.var,
main.vars=main.vars, int.vars=int.vars, strata.var=strata.var, ProbG1.var=ProbG1.var,
op=op), silent=TRUE)
errorFlag <- ("try-error" %in% class(temp))
if (!errorFlag) {
if (firstFlag) {
writeOut(fid, c(snp.name, majMin, maf, temp[-1]))
} else {
writeOut(fid, c("SNP", "Alleles", "MAF", names(temp[-1])))
writeOut(fid, c(snp.name, majMin, maf, temp[-1]))
firstFlag <- 1
errorVec <- temp[-1]
errorVec[] <- "NA"
}
} else {
if (firstFlag) writeOut(fid, c(snp.name, majMin, maf, errorVec))
}
} # END: while(1)
close(fid)
op$out.file
} # END: scan.UML_CML
# Function for meta-anlaysis
meta.fixed <- function(study.list) {
# study.list List of named lists with names beta and cov
nstudy <- length(study.list)
sigma.inv.list <- list()
sum.sigma.inv <- 0
sum.sigma.inv.beta <- 0
# Compute sigma inverse for each study
for (i in 1:nstudy) {
temp <- study.list[[i]]
beta <- temp$beta
sigma <- temp$cov
sigma.inv <- solve(sigma)
sum.sigma.inv <- sum.sigma.inv + sigma.inv
dim(beta) <- c(length(beta), 1)
sum.sigma.inv.beta <- sum.sigma.inv.beta + sum.sigma.inv %*% beta
sigma.inv.list[[i]] <- sigma.inv
}
meta.cov <- solve(sum.sigma.inv)
meta.beta <- meta.cov %*% sum.sigma.inv.beta
list(meta.beta=meta.beta, meta.cov=meta.cov, meta.cov.inv=sum.sigma.inv,
inv.list=sigma.inv.list)
} # END: meta.fixed
# Function to compute A matrix for EB meta
get.EB.A_matrix <- function(parm1, cov1, parm2) {
psi <- parm1 - parm2
dim(psi) <- NULL
psi2 <- psi*psi
phi <- diag(cov1)
temp <- psi2/(psi2 + phi)
A <- diag(temp)
A
} # END: get.EB.A_matrix
meta.EB <- function(beta1, beta2, sumInv1.inv, sumInv2.inv, inv1.list, inv2.list,
cov12.list) {
# Covariance matrix for beta1 is sumInv1.inv
A <- get.EB.A_matrix(beta1, sumInv1.inv, beta2)
n <- length(beta1)
I <- matrix(0, nrow=n, ncol=n)
diag(I) <- 1
beta.EB <- (A %*% beta1) + (I - A) %*% beta2
sum <- 0
n.study <- length(inv1.list)
for (i in 1:n.study) {
sum <- sum + inv1.list[[i]] %*% cov12.list[[i]] %*% inv2.list[[i]]
}
cov.EB <- sumInv1.inv %*% sum %*% sumInv2.inv
list(beta=beta.EB, cov=cov.EB)
} # END: meta.EB
# Function to get the parms/cov for a study given a named vec
getBetaObject <- function(vec, vars) {
# vec Named vector
# vars Names to extract (must be ordered G, E, GE)
ret <- vec[vars]
ret
} # END: getBetaObject
# Function to get the cov for a study given a named vec
getCovObject <- function(vec, vars) {
# vec Named vector
# vars Names to extract (must be ordered, upper triangle of matrix)
len <- length(vars)
n <- (-1 + sqrt(8*len + 1))/2
ret <- matrix(data=NA, nrow=n, ncol=n)
a <- 1
for (i in 1:n) {
b <- a + n - i
ret[i, i:n] <- vec[vars[a:b]]
ret[i:n, i] <- ret[i, i:n]
a <- b + 1
}
ret
} # END: getCovObject
# Function to get the UML-CML cov for a study given a named vec
getCov12Object <- function(vec, vars) {
# vec Named vector
# vars Names to extract (must be ordered, rows of matrix)
n <- sqrt(length(vars))
ret <- matrix(data=vec[vars], byrow=TRUE, nrow=n, ncol=n)
ret
} # END: getCov12Object
# Function to perfrm meta-analysis for each study
meta.GxE <- function(study.list, op=NULL) {
# study.list List of sublists containing name of file, study name, etc
# op List with names:
# out.file Output file
# snps List of snps to include
# print 0 or 1
##############################################################
# Function to get the number of E parms
##############################################################
get.n.parms <- function(vec, which) {
# vec character vector of column names
# which "E" (for main effects) or "G" (for interactions)
# Remove UML.G.BETA and CML.G.BETA
temp <- !(vec %in% c("UML.G.BETA", "CML.G.BETA"))
vec <- vec[temp]
ids <- grep("BETA", vec, fixed=TRUE)
vec <- vec[ids]
vec <- gsub("UML.", "", vec, fixed=TRUE)
vec <- gsub("CML.", "", vec, fixed=TRUE)
vec <- gsub(".BETA", "", vec, fixed=TRUE)
temp <- substr(vec, 1, 1) == which
vec <- vec[temp]
vec <- unique(vec)
n <- length(vec)
n
} # END: get.n.parms
###############################################################
# Function to get the parm names that would need to be flipped from
# a different reference cat
get.flip.E <- function(vec) {
# Remove UML.G.BETA and CML.G.BETA
temp <- !(vec %in% c("UML.G.BETA", "CML.G.BETA"))
vec <- vec[temp]
ids <- grep("BETA", vec, fixed=TRUE)
vec <- vec[ids]
vec2 <- vec
vec2 <- gsub("UML.", "", vec2, fixed=TRUE)
vec2 <- gsub("CML.", "", vec2, fixed=TRUE)
temp <- substr(vec2, 1, 1) %in% c("E","G")
ret <- vec[temp]
ret
} # END: get.flip.E
###############################################################
# Function to determine which column names would be flipped from
# a different risk allele
get.flip.cnames <- function(vec) {
ids <- grep("BETA", vec, fixed=TRUE)
vec <- vec[ids]
vec2 <- vec
vec2 <- gsub("UML.", "", vec2, fixed=TRUE)
vec2 <- gsub("CML.", "", vec2, fixed=TRUE)
temp <- substr(vec2, 1, 1) == "G"
ret <- vec[temp]
ret
} # END: get.flip.cnames
#################################################################
# Function to flip an allele
flip.allele <- function(a) {
if (a == "A") {
return("T")
} else if (a == "T") {
return("A")
} else if (a == "C") {
return("G")
} else if (a == "G") {
return("C")
} else {
return(NULL)
}
} # END: flip.allele
#################################################################
# Function to determine if G parms need to be flipped
get.flip.status <- function(risk.base, other.base, MAF.base,
risk, other, MAF) {
message <- ""
messageFlag <- 0
flip <- 0
error <- 0
# See if the set of alleles match
set.base <- c(risk.base, other.base)
set <- c(risk, other)
flag <- all(set %in% set.base)
if (!flag) {
# Try flipping set of alleles
set[1] <- flip.allele(set[1])
set[2] <- flip.allele(set[2])
flag <- all(set %in% set.base)
}
if (!flag) {
message <- "Alleles do not match"
messageFlag <- 1
return(flip=0, error=1, message=message)
}
# See if risk allele matches
if (risk != set[1]) flip <- 1
list(flip=flip, error=error, message=message)
} # END: get.flip.status
###################################################################
op <- default.list(op, c("out.file", "print"), list("ERROR", 0),
error=c(1, 0))
nstudy <- length(study.list)
snps <- op[["snps", exact=TRUE]]
snpFlag <- !is.null(snps)
snps.all <- NULL
cnames <- NULL
outfile <- op[["out.file", exact=TRUE]]
outFlag <- !is.null(outfile)
if (snpFlag) snps <- removeWhiteSpace(snps)
# Check each study list
vec <- scan(study.list[[1]]$file, what="character", nlines=1, quiet=TRUE)
for (i in 1:nstudy) {
study.list[[i]] <- check.file.list(study.list[[i]], op=list(vars=vec))
temp <- study.list[[i]][["flip.betas", exact=TRUE]]
if (is.null(temp)) study.list[[i]]$flip.betas <- 0
}
# Read in the data
for (i in 1:nstudy) {
tlist <- study.list[[i]]
x <- loadData.table(tlist)
x[, "SNP"] <- removeWhiteSpace(x[, "SNP"])
if (snpFlag) {
temp <- x[, "SNP"] %in% snps
x <- removeOrKeepRows(x, temp)
}
nx <- nrow(x)
str <- paste("For study ", i, ", ", nx, " SNPs will be included.", sep="")
print(str)
if (!nx) {
str <- paste("Removing study ", i, sep="")
print(str)
study.list[[i]] <- NULL
next
}
rownames(x) <- x[, "SNP"]
tlist$data <- x
study.list[[i]] <- tlist
snps.all <- unique(c(snps.all, x[, "SNP"]))
# Check the names of the columns
if (i == 1) {
cnames <- colnames(x)
} else {
if (!all(cnames == colnames(x))) {
str <- paste("ERROR with column names for study ", i, sep="")
stop(str)
}
}
}
# Get the number of main effect parms and interaction parms
N.E.parms <- get.n.parms(cnames, "E")
N.GE.parms <- get.n.parms(cnames, "G")
# Get the column names for the vector of main effect estimates
UML.beta.cnames <- outNames.me("UML", N.E.parms, N.GE.parms)
CML.beta.cnames <- outNames.me("CML", N.E.parms, N.GE.parms)
UML.cov.cnames <- outNames.cov("UML", N.E.parms, N.GE.parms)
CML.cov.cnames <- outNames.cov("CML", N.E.parms, N.GE.parms)
UML.CML.cov.cnames <- outNames.cov2(N.E.parms, N.GE.parms)
# Get the betas to flip from different reference cat
cat.flip.cnames <- get.flip.E(cnames)
# Get the beta column names to flip for different risk allele
allele.flip.cnames <- get.flip.cnames(cnames)
outnames <- c("SNP", "Risk.allele", "Other.allele", "P.omnibus",
"N.study", "Studies", "Message")
outvec <- character(length(outnames))
names(outvec) <- outnames
# Open the output file
if (outFlag) {
fid <- file(outfile, "w")
writeOut(fid, outnames)
}
# Loop over each snp
nsnps <- length(snps.all)
for (i in 1:nsnps) {
outvec[] <- "NA"
UML.list <- list()
CML.list <- list()
index <- 0
snp <- snps.all[i]
message <- ""
outvec[1:3] <- c(snp, ".", ".")
# For UML-CML cov matrix
UML.CML.list <- list()
# Keep track of alleles and MAF
risk.vec <- NULL
other.vec <- NULL
MAF.vec <- NULL
study.vec <- NULL
# Get the vector of estimates for each study
for (j in 1:nstudy) {
x <- study.list[[j]]$data
if (!(snp %in% rownames(x))) next
vec <- x[snp, ]
names(vec) <- cnames
alleles <- vec["Alleles"]
other <- substr(alleles, 1, 1)
risk <- substr(alleles, 2, 2)
MAF <- as.numeric(vec["MAF"])
vec <- as.numeric(vec[-(1:3)])
if (any(!is.finite(vec))) next
names(vec) <- cnames[-(1:3)]
# Flip betas if needed
if (study.list[[j]]$flip.betas) vec[cat.flip.cnames] <- -vec[cat.flip.cnames]
# Get the UML and CML estimates
CML.beta <- getBetaObject(vec, CML.beta.cnames)
CML.cov <- getCovObject(vec, CML.cov.cnames)
UML.beta <- getBetaObject(vec, UML.beta.cnames)
UML.cov <- getCovObject(vec, UML.cov.cnames)
# For UML-CML cov matrix
UML.CML.cov <- getCov12Object(vec, UML.CML.cov.cnames)
if (any(!is.finite(UML.CML.cov))) next
# Check the alleles
if (index == 0) {
risk.base <- risk
other.base <- other
outvec[2:3] <- c(risk, other)
MAF.base <- MAF
} else {
# Compare to base values
ret <- get.flip.status(risk.base, other.base, MAF.base, risk, other, MAF)
if (ret$error) next
if (ret$flip) {
# Flip betas
UML.beta[allele.flip.cnames] <- -UML.beta[allele.flip.cnames]
CML.beta[allele.flip.cnames] <- -CML.beta[allele.flip.cnames]
}
}
index <- index + 1
#alleles.vec <- c(alleles.vec, alleles)
#MAF.vec <- c(MAF.vec, MAF)
study.vec <- c(study.vec, j)
# Add parms to the lists
tlist <- list(beta=UML.beta, cov=UML.cov)
UML.list[[index]] <- tlist
tlist <- list(beta=CML.beta, cov=CML.cov)
CML.list[[index]] <- tlist
tlist <- list(beta=CML.beta, cov=CML.cov)
UML.CML.list[[index]] <- UML.CML.cov
}
outvec["N.study"] <- index
if (!index) {
outvec["Message"] <- "No studies"
if (outFlag) writeOut(fid, outvec)
next
}
outvec["Studies"] <- paste(study.vec, collapse=",", sep="")
# Perform the fixed effects meta-anlaysis for UML and CML estimates
# Return list: list(meta.beta=meta.beta, meta.cov=meta.cov, meta.cov.inv=sum.sigma.inv,
# inv.list=sigma.inv)
UML.meta <- try(meta.fixed(UML.list), silent=TRUE)
if (checkTryError(UML.meta, conv=0)) {
outvec["Message"] <- "Error in UML meta-analysis"
if (outFlag) writeOut(fid, outvec)
next
}
CML.meta <- try(meta.fixed(CML.list), silent=TRUE)
if (checkTryError(CML.meta, conv=0)) {
outvec["Message"] <- "Error in CML meta-analysis"
if (outFlag) writeOut(fid, outvec)
next
}
# Now perform EB
# meta.EB <- function(beta1, beta2, sumInv1.inv, sumInv2.inv, inv1.list, inv2.list,
# cov12.list)
# Return list: list(beta=beta.EB, cov=cov.EB)
ret <- meta.EB(UML.meta$meta.beta, CML.meta$meta.beta, UML.meta$meta.cov,
CML.meta$meta.cov, UML.meta$inv.list,
CML.meta$inv.list, UML.CML.list)
temp <- waldTest.main(ret$beta, ret$cov, 1:length(UML.beta))
outvec["P.omnibus"] <- temp$pvalue
outvec["Message"] <- "OK"
# Write output
if (outFlag) writeOut(fid, outvec)
}
if (outFlag) close(fid)
outvec
} # END: meta.GxE
Try the CGEN package in your browser
Any scripts or data that you put into this service are public.
CGEN documentation built on April 28, 2020, 8:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions meta.GxE getCov12Object getCovObject getBetaObject meta.EB get.EB.A_matrix meta.fixed scan.UML_CML extractEst extract_UML.CML outNames.cov2 outNames.cov outNames.me writeOut UML_CML_GxE_parms apply_zero.vars logistic.dsgnMat logistic.vnames snp.main snp.logistic snp.scan.logistic sMatrix.logistic
Documented in snp.logistic
# History Apr 21 2008 Initial coding
# May 09 2008 Use returned hessian to compute cov
# May 29 2008 Redo the log-likelihood calculation
# May 31 2008 Move tests to GbyE function.
# Jun 02 2008 Add gradient function for the optimizer
# Add code for new output file.
# Jun 03 2008 Add test2.vars for a user specified Wald test.
# Jun 04 2008 Add tests.init (standard logistic reg.)
# Jun 06 2008 Add parameter estimates and standard errors
# to out.file.
# Make vnames a seperate function
# Jun 07 2008 Combine GbyE.1snp and GbyE.irls
# Jun 09 2008 Change input argument names, function names.
# Remove strataMat from input.
# Jun 13 2008 Catch errors when calling logistic.SNP
# Jun 16 2008 Rename getLoglike to getLoglike.glm
# Jun 17 2008 Rename strata vars to Allele_freq
# Put strata parms in a separate object
# Jun 18 2008 Update for change in getWaldTest
# Jun 20 2008 Change X.vars to X.main, V.vars to X.int
# Jun 30 2008 Remove call to setUpSummary
# Jul 02 2008 Add code for streamed input
# Jul 09 2008 Change in getEB function:
# cmat <- cbind(diag(temp), diag(1-temp))
# Jul 11 2008 Rename the functions to snp.logistic and
# snp.scan.logistic
# Jul 11 2008 Add output for snp and subject counts
# Jul 17 2008 Add code for factor snps
# Jul 18 2008 Add test.inter
# Jul 18 2008 Add code for saving individual tests and p-values
# Jul 21 2008 Add option to output corr parms
# Jul 22 2008 Change code to call getTest function
# Jul 26 2008 Add tests and tests.1df options
# The tests option will replace the test2.vars option
# Aug 07 2008 Check for constant variables in getDsgnMat
# Aug 11 2008 Define classes for return objects
# Aug 21 2008 Add option for joint and stratified effects
# Aug 25 2008 Rename snp.logistic to snp.main
# Add snp.logistic function
# Aug 29 2008 Update for change in getWaldTest
# Oct 08 2008 Allow for different family if only.UML = 1
# Oct 10 2008 Allow for no recoding of genotypes when stream = 1
# Oct 21 2008 Make more efficient for only.UML = 1
# Oct 29 2008 Change the line in snp.scan.logistic from
# snp <- as.integer(snp) to
# snp <- as.numeric(snp)
# Dec 10 2008 Allow the user to specify formulas in
# snp.logistic.
# Jan 30 2009 Add option to output UML base model
# Mar 05 2009 Get MAF from getData.1 function (stream = 0)
# Mar 05 2009 Update for stream = 1
# Mar 16 2009 Do not set cc.var to response.var by default
# Jun 16 2009 Print out response levels only for binomial family
# Jul 10 2009 Check if design.V0 is NULL in getInter.Vars()
# Jul 22 2009 Move temp.list inside the options list
# Aug 03 2009 Change in outputRow to match the change in
# effects.init return list
# Sep 17 2009 Add in output option in getDelimiter in snp.scan.logistic
# Change as.integer to as.numeric
# Add errorCheck option in snp.main
# Oct 01 2009 Check variable names in snp.logistic
# Oct 16 2009 Change in recode.geno
# Oct 16 2009 Add major/minor allele to out file in snp.scan.logistic
# Oct 20 2009 Fix bug in getting response0 in snp.scan.logistic
# Sometimes it is character
# Oct 28 2009 Make sure design matrices from logistic.dsgnMat
# are numeric
# Dec 17 2009 Remove IRLS algorithm
# Add code for genetic.model = 3
# Dec 23 2009 Generalize the stratification
# Jan 04 2010 Update snp.logistic for missing values
# Feb 01 2010 Wrap callOptim in try function in snp.main
# Feb 04 2010 In snp.logistic, check for 1 strata.var that is constant
# Feb 12 2010 Use chol and chol2inv instead of solve.
# Add option for C code.
# Feb 17 2010 Check CML variances from returned C code
# Feb 17 2010 Make getInit function in snp.main more efficient
# Feb 26 2010 Use the snp variable name in snp.logistic
# Mar 04 2010 Use normVarNames function in snp.logistic
# Mar 14 2010 Add option for using (2) only controls to determine
# the major/minor allele or (1) all subjects
# The default is 1
# Mar 18 2010 Add package="CGEN" option
# Mar 29 2010 Add return values to snp.logistic
# Mar 30 2010 Add option for fixing parameters
# Apr 02 2010 Add option zero.vars to snp.logistic
# Apr 05 2010 Change call to C code
# Apr 14 2010 Change default reltol to 1e-8
# Apr 22 2010 Fix bug in getInit with adding snp back when zeroSNP = 1
# Apr 26 2010 Change code for when there are NAs with UML parms
# Nov 04 2011 Check if snp.var is a main effect
# Dec 23 2011 Add function in snp.main to check CML initial estimates
# Allow for initial estimates to be passed in
# Jul 18 2012 Remove errorCheck option in snp.main. Was causing a problem
# for different genetic models.
# Feb 02 2013 Add warning message for continuous variables
# Mar 26 2013 Add code for imputed SNPs
# Mar 29 2013 Add new functions for outputing UML-CML estimates
# Apr 03 2013 Add code for meta-analysis
# Apr 04 2013 Generalize code for UML-CML estimates
# Apr 05 2013 Modify snp.scan.logistic for imputed snps
# Sep 01 2013 Allow snp to be a 3 column matrix in snp.main
# Oct 22 2013 Return full CML covariance matrix from snp.main
# Oct 25 2013 Add function to return strat matrix
######################################################################
# TO DO:
# Modify the different tests for the different genetic models
# Update test2 code for factors
# Check for consistent input arguments (out.est, only.UML,...)
######################################################################
sMatrix.logistic <- function(data, strata.var, facVars) {
sflag <- !is.null(strata.var)
sform <- "formula" %in% class(strata.var)
nobs <- nrow(data)
svarCat1 <- FALSE
# Check for constant strata variable
if ((sflag) && (!sform) && (length(strata.var) == 1)) {
svarCat1.v <- data[, strata.var]
svarCat1.s <- sort(unique(svarCat1.v))
svarCat1.n <- length(svarCat1.s)
if (svarCat1.n == 1) sflag <- FALSE
if (sflag) {
if ((is.character(svarCat1.v)) || (is.factor(svarCat1.v))) svarCat1 <- TRUE
}
}
if (!sflag) {
design.S0 <- matrix(data=1, nrow=nobs, ncol=1)
s.1catVar <- 1
} else {
if (svarCat1) {
design.S0 <- matrix(data=0, nrow=nobs, ncol=svarCat1.n)
s.1catVar <- 1
for (i in 1:svarCat1.n) design.S0[(svarCat1.v %in% svarCat1.s[i]), i] <- 1
} else {
design.S0 <- logistic.dsgnMat(data, strata.var, facVars, removeInt=0)$designMatrix
s.1catVar <- 0
}
}
list(design.S0=design.S0, s.1catVar=s.1catVar)
} # END: sMatrix.logistic
# Function to perform SNP by environment interaction analysis.
snp.scan.logistic <- function(snp.list, pheno.list, op=NULL) {
# INPUT:
# snp.list (See snp.list.wordpad)
#################################################################
# pheno.list (See pheno.list.wordpad)
# This list must include the name "response.var",
# and can also include the names "strata.var",
# "main.vars", "int.vars" (see below)
# response.var Name of the response variable. This variable
# must be coded as 0 and 1.
# No default
# strata.var Stratification variable name(s) or a formula.
# The default is NULL so that all observations
# belong to the same strata.
# main.vars Vector of variables names in the file pheno.list$file
# that will be included in the model as main effects.
# main.vars can also be a formula.
# The default is NULL.
# int.vars Vector of variable names in the file pheno.list$file
# that will interact with each snp in the snp data.
# int.vars can also be a formula.
# The default is NULL.
# cc.var The name of the case-control variable. This variable must
# be coded as 0-1. The variable is used to determine
# the MAF based on controls.
# The default is NULL.
#################################################################
# op List with the following names.
# genetic.model 0-2
# 0: trend
# 1: dominant
# 2: recessive
# 3: general
# The default is 0.
# reltol Stopping tolerance
# The default is 1e-8
# maxiter Maximum number of iterations
# The default is 100
# optimizer One of : "BFGS", "Nelder-Mead", "BFGS", "CG",
# "L-BFGS-B", "SANN".
# The default is "BFGS"
# print 0 or 1 to print each list of output
# The default is 0
# test.omnibus 0 or 1 for the test
# involving the snp and interaction terms.
# The default is 1
# test.main 0 or 1 for computing
# the test for the main effect of the snp.
# The default is 1.
# test.inter 0 or 1 for computing
# the test for the interaction terms.
# The default is 0.
# tests List of character vector of variable names that
# will be used in Wald tests.
# The default is NULL.
# tests.1df Character vector of variable names to use for 1 df
# Wald tests.
# The default is NULL.
# tests.UML NULL or a list with at least one of the following
# names: "omnibus", "main", or "inter". These names
# can be set to one of the following values:
# "WALD", "LR", or "SCORE".
# Example: test.UML=list(main="WALD", omnibus="LR") will
# compute a wald test for the main effect of SNP and
# a likelihood ratio test for the SNP and interactions.
# !!!! Only for family=binomial !!!!!!
# The default value for each name is "WALD"
# The default is NULL.
# tests.CML 0 or 1 for computing the WALD tests using the
# conditional ML estimates.
# The default is 1.
# tests.EB 0 or 1 for computing the WALD tests using the
# empirical bayes estimates
# The default is 0.
# geno.counts 0 or 1 to add the genotype counts to out.file
# The default is 0.
# subject.counts 0 or 1 to add the number of cases and controls to
# out.file.
# The default is 0.
# allele.cc 1 or 2 to use all subjects
###############################################################
# effects List for joint/stratified effects
# Names in the list must be "var", "type",
# "var.levels", "snp.levels", "var.base"
# The default is NULL.
# var Variable name to compute the effects with the SNP
# No default
# type 1 or 2 or c(1, 2) , 1 = joint, 2 = stratified
# The default is 1.
# var.levels (Only for continuous var) Numeric vector of the
# levels to be used in the calculation.
# The default is 0.
# var.base (Only for continuous var) Baseline level.
# The default is 0
# snp.levels One of the following: 0, 1, 2, c(0,1), c(0,2),
# c(1,2) or c(0, 1, 2)
# The default is 1.
# method Character vector containing any of the following:
# "UML", "CML", "EB"
# The default is c("UML", "CML", "EB")
###############################################################
# out.est List with the names "parms" and "method".
# Use this option if you want to save parameter
# estimates and standard errors to out.file.
# parms 1-3 or a character vector of the parameters to
# save statistics on.
# 1: Main effect is saved.
# 2: Main effect and interactions are saved.
# 3: All parameters are saved.
# method Character vector containing any of the following:
# "UML", "CML", "EB"
# For example, setting parms=3 and method=c("CML", "EB")
# will write all the parameter estimates and
# standard errors for the methods CML and EB.
# The default value of out.est is NULL.
# what Character vector containing any of the following:
# "beta", "se", "test", "pvalue"
# corr.parms List of character vectors of length 2 giving the
# pairs of variables for correlations.
###############################################################
# out.file NULL or file name to save summary information for
# each snp. The output will at least contain the columns
# "SNP" and "MAF". MAF is the minor allele frequency
# form the controls. Additional columns in this file
# are based on the values of test.omnibus, test.main,
# test2.vars, tests.UML, tests.CML, tests.EB, and
# out.est.
# The default is NULL.
# out.dir NULL or the output directory to store the output
# lists for each SNP. A seperate file will be created
# for each snp in the snp data set.
# The file names will be the out_<snp>.rda
# The object names are called "ret".
# The default is NULL.
# base.outfile NULL or path to base model rda file
# The default is NULL.
#####################################################################
# temp.list See temp.list.doc
#####################################################################
# RETURN:
# The list returned is from the last analysis performed.
# Function to call different test functions
getTest <- function(fit, vars, test, which, model) {
# fit List of parms, cov, and loglike
# vars Variables to test
# test "WALD", "LR", or "SCORE"
# which 1-4 1=omnibus, 2=main, 3=inter, 4=test2
# model "UML", "CML", or "EB"
if (test == "WALD") {
# Set up a list for calling getWaldTest
temp <- list(parms=fit$parms, cov=fit$cov)
return(getWaldTest(temp, vars))
}
# Determine if test is valid
if (test == "LR") {
# LR is not for EB
if (model == "EB") return(list(test=NA, df=NA, pvalue=NA))
if ((model == "CML") && (which %in% c(1, 2))) return(list(test=NA, df=NA, pvalue=NA))
} else if (test == "SCORE") {
if ((model == "CML") && (which %in% c(1, 2))) return(list(test=NA, df=NA, pvalue=NA))
}
# Determine if interactions are in the model
if (which == 2) {
X.int <- design.V
xflag <- 1
} else {
X.int <- NULL
xflag <- 0
}
# Factor snp if genetic.model = 3
if (genetic.model == 3) snp <- factor(snp)
# Determine if SNP is in the model
if (which == 3) {
int.vec <- snp
flag <- 1
} else {
int.vec <- NULL
flag <- 0
}
# Base model
if (model == "UML") {
fit0 <- callGLM(response, X.main=design.X, X.int=X.int, int.vec=snp,
family="binomial", prefix="SNP_", retX=TRUE,
retY=TRUE, inc.int.vec=flag)
} else {
fit0 <- snp.main(response, snp, X.main=design.X, X.int=design.V,
X.strata=design.S, op=op)
fit0 <- fit0[[model]]
# Add info for score tests
if (test == "SCORE") {
# Add x, y, linear.predictors
temp <- getXBeta(cbind(design.X, int.vec, X.int), fit0$parms)
fit0$linear.predictors <- temp$XBeta
fit0$x <- temp$X
fit0$y <- response
} else {
# LR test, add the rank
fit0$rank <- ncol(fit0$cov)
}
}
if (test == "LR") {
# Likelihood ratio test. Add the rank
fit$rank <- ncol(fit$cov)
return(likelihoodRatio(fit, fit0))
} else {
# Score test
# Get the matrix to test
temp <- NULL
if ((genetic.model == 3) && (!xflag || !flag)) {
snp <- createDummy(snp)$data
}
if (!xflag) temp <- addInterVars(NULL, snp, design.V)$data
if (!flag) temp <- cbind(temp, snp)
return(score.logReg(fit0, temp))
}
} # END: getTest
# Function to add tests to the return list
addToList <- function(ret, which) {
# ret Return list
# which "UML", "CML", "EB"
if (which == "UML") {
# tests.UML.vec is actually a list
test <- tests.UML.vec
} else {
test <- rep("WALD", times=3)
}
field <- paste(which, c(".omnibus", ".main", ".inter"), sep="")
temp <- ret[[which]]
if (!is.null(temp)) {
if (omniFlag) ret[[field[1]]] <-
getTest(temp, omni.vars[[vListIndex]], test[1], 1, which)
if (mainFlag) ret[[field[2]]] <-
getTest(temp, main.vars[[vListIndex]], test[2], 2, which)
if (interFlag) ret[[field[3]]] <-
getTest(temp, inter.vars[[vListIndex]], test[3], 3, which)
# tests will start from field 4
if (test2Flag) {
field2 <- paste(which, ".test", 1:n.tests, sep="")
for (i in 1:n.tests) {
ret[[field2[i]]] <- getTest(temp, test2.vars[[i]], "WALD", 4, which)
}
}
} else {
if (omniFlag) ret[[field[1]]] <- list(test=NA, pvalue=NA, df=NA)
if (mainFlag) ret[[field[2]]] <- list(test=NA, pvalue=NA, df=NA)
if (interFlag) ret[[field[3]]] <- list(test=NA, pvalue=NA, df=NA)
if (test2Flag) {
temp <- list(test=NA, pvalue=NA, df=NA)
for (i in 1:n.tests) ret[[field2[i]]] <- temp
}
}
ret
} # END: addTest
# Function to tests to the output vector
addToVec <- function(outVec, which) {
# outVec Output vector
# which "omnibus", "main", "inter", or "testi"
tname <- paste(which, ".test", sep="")
pname <- paste(which, ".pvalue", sep="")
dname <- paste(which, ".df", sep="")
if (tests.UML) {
t <- paste("UML.", tname, sep="")
p <- paste("UML.", pname, sep="")
d <- paste("UML.", dname, sep="")
r <- paste("UML.", which, sep="")
outVec[t] <- ret[[r]]$test
outVec[p] <- ret[[r]]$pvalue
outVec[d] <- ret[[r]]$df
}
if (tests.CML) {
t <- paste("CML.", tname, sep="")
p <- paste("CML.", pname, sep="")
d <- paste("CML.", dname, sep="")
r <- paste("CML.", which, sep="")
outVec[t] <- ret[[r]]$test
outVec[p] <- ret[[r]]$pvalue
outVec[d] <- ret[[r]]$df
}
if (tests.EB) {
t <- paste("EB.", tname, sep="")
p <- paste("EB.", pname, sep="")
d <- paste("EB.", dname, sep="")
r <- paste("EB.", which, sep="")
outVec[t] <- ret[[r]]$test
outVec[p] <- ret[[r]]$pvalue
outVec[d] <- ret[[r]]$df
}
outVec
} # END: addToVec
# Function for the number of cases and controls
getCCcounts <- function(notMiss) {
# Note: The cntrl vector equals 1 for controls.
# So the returned vector from table(), will have the counts
# for the cases and then the controls, which is the order in
# the output file.
if (is.null(notMiss)) return(subj.cnts0)
temp <- cntrl[notMiss]
nn <- sum(temp)
return(c(length(temp)-nn, nn))
} # END: getCCcounts
# Function to get variables for main test
getMain.vars <- function() {
v2 <- NULL
v1 <- getVarNames.snp(prefix=op$snpName, genetic.model=genetic.model)
if (genetic.model == 3) {
v2 <- getVarNames.snp(prefix=op$snpName, genetic.model=0)
}
list(v1, v2)
} # END: getMain.vars
# Function to get variables for omnibus test
getOmni.vars <- function() {
# Combine main effect and interaction vars
mvars <- getMain.vars()
temp <- getInter.vars()
v1 <- c(mvars[[1]], temp[[1]])
v2 <- c(mvars[[2]], temp[[2]])
list(v1, v2)
} # END: getMain.vars
# Function to get the interaction vars
getInter.vars <- function() {
if (is.null(design.V0)) return(NULL)
v2 <- NULL
v1 <- getVarNames.int(design.V0, prefix=op$snpName, genetic.model=genetic.model, sep=":")
if (genetic.model == 3) {
v2 <- getVarNames.int(design.V0, prefix=op$snpName, genetic.model=0, sep=":")
}
list(v1, v2)
} # END: getInter.vars
# Function to return the test2 vars
getTest2.vars <- function() {
tests <- getListName(op, "tests")
ret <- list()
# Loop over each set of variables
for (i in 1:length(tests)) {
testi <- tests[[i]]
# Determine if any of the test2 variables are factors
temp <- 0
if (facFlag) temp <- testi %in% facVars
if (!any(temp)) {
test2.vars <- testi
} else {
# First add the variables that are not factors
test2.vars <- testi[as.logical(1-temp)]
# Now add the dummy variables for the factors
temp <- testi[temp]
for (temp2 in temp) {
# Get the dummy variable names in one of the lists
temp3 <- getListName(X.newVars, temp2)
if (!is.null(temp3)) {
test2.vars <- c(test2.vars, temp3)
} else {
temp3 <- getListName(V.newVars, temp2)
test2.vars <- c(test2.vars, temp3)
}
}
}
# Check the variable names
temp <- check.vec.char(test2.vars, "tests", len=NULL,
checkList=log.vnames$all)
if (temp) stop("ERROR with op$tests")
# Add to the list
ret[[i]] <- test2.vars
}
ret
} # END: getTest2.vars
# Function to return the index to use to get the test variable names
getVListIndex <- function(n) {
if (genetic.model != 3) {
return(1)
} else {
if (n == 3) {
return(1)
} else {
return(2)
}
}
} # END: getVListIndex
# Function to in initialize the output file
initOutfile <- function() {
# Create an output vector
vnames <- NULL
vec <- c(".test", ".pvalue", ".df")
if (omniFlag) {
temp <- paste("omnibus", vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
if (mainFlag) {
temp <- paste("main", vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
if (interFlag) {
temp <- paste("inter", vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
if (test2Flag) {
for (i in 1:n.tests) {
temp <- paste("test", i, vec, sep="")
if (tests.CML) vnames <- c(vnames, paste("CML.", temp, sep=""))
if (tests.UML) vnames <- c(vnames, paste("UML.", temp, sep=""))
if (tests.EB ) vnames <- c(vnames, paste("EB.", temp, sep=""))
}
}
est.se <- NULL
est.test <- NULL
est.pval <- NULL
est.corr <- NULL
if (out.est.flag) {
if ((out.est$parms == 2) && (!Vflag)) out.est$parms <- 1
# Vector to hold parm names
what <- out.est$what
temp <- NULL
if (out.est.beta) temp <- est.p[[1]]
if (out.est.se) {
est.se <- list()
for (i in 1:length(est.parms)) {
est.se[[i]] <- paste(est.parms[[i]], ".se", sep="")
}
temp <- c(temp, est.se[[1]])
}
if (out.est.test) {
est.test <- list()
for (i in 1:length(est.parms)) {
est.test[[i]] <- paste(est.parms[[i]], ".test", sep="")
}
temp <- c(temp, est.test[[1]])
}
if (out.est.pval) {
est.pval <- list()
for (i in 1:length(est.parms)) {
est.pval[[i]] <- paste(est.parms[[i]], ".pvalue", sep="")
}
temp <- c(temp, est.pval[[1]])
}
if (out.est.corr) {
est.corr <- NULL
corrs <- getListName(out.est, "corr.parms")
for (i in 1:length(corrs)) {
var <- paste(corrs[[i]], collapse="_", sep="")
est.corr <- c(est.corr, var)
}
temp <- c(temp, est.corr)
}
for (method in out.est$method) {
temp2 <- paste(method, ".", temp, sep="")
vnames <- c(vnames, temp2)
}
}
# Effects
if (effectsFlag) {
names <- effects$out.names
nr <- nrow(names)
for (method in effects$method) {
for (type in effects$type) {
for (i in 1:nr) {
temp <- paste(method, ".eff", type, ".",
names[i,], sep="")
vnames <- c(vnames, temp)
}
# Standard errors
for (i in 1:nr) {
temp <- paste(method, ".eff", type, ".",
names[i,], ".se", sep="")
vnames <- c(vnames, temp)
}
}
}
}
outVec <- double(length(vnames))
names(outVec) <- vnames
# geno and subject counts
if (subj.counts) vnames <- c("N.cases", "N.controls", vnames)
if (geno.counts) vnames <- c("Hom.common", "Heter", "Hom.uncommon", vnames)
vnames <- c("SNP", "Alleles", "MAF", vnames)
fid <- writeVecToFile(vnames, op$out.file, colnames=NULL, type=3,
close=0, sep="\t")
list(fid=fid, outVec=outVec, est.se=est.se,
est.test=est.test, est.pval=est.pval, est.corr=est.corr)
} # END: initOutfile
# Function to create a design matrix
getDsgnMat <- function(vars) {
design <- removeOrKeepCols(phenoData0, vars, which=1)
newVars <- NULL
if (facFlag) {
temp <- vars %in% facVars
if (any(temp)) {
temp <- vars[temp]
temp <- createDummy(design, vars=temp)
design <- temp$data
newVars <- temp$newVars
}
}
design <- as.matrix(design)
# Check for constant variables
design <- checkForConstantVar(design, msg=1)$data
list(designMatrix=design, newVars=newVars)
} # END: getDsgnMat
# Function to output a row to the output file
outputRow <- function() {
cat(snp.name, majMin, maf, "", file=fid, sep="\t")
if (geno.counts) cat(genoCounts, "", file=fid, sep="\t")
if (subj.counts) cat(getCCcounts(subj.notMiss), "", file=fid, sep="\t")
if (omniFlag) outVec <- addToVec(outVec, "omnibus")
if (mainFlag) outVec <- addToVec(outVec, "main")
if (interFlag) outVec <- addToVec(outVec, "inter")
if (test2Flag) {
for (i in 1:n.tests) {
outVec <- addToVec(outVec, paste("test", i, sep=""))
}
}
if (out.est.flag) {
parms <- est.p[[vListIndex]]
for (method in out.est$method) {
temp <- getListName(ret, method)
if (out.est.beta) temp1 <- paste(method, ".", est.parms[[vListIndex]], sep="")
if (out.est.se) temp2 <- paste(method, ".", est.se[[vListIndex]], sep="")
if (out.est.test) temp3 <- paste(method, ".", est.test[[vListIndex]], sep="")
if (out.est.pval) temp4 <- paste(method, ".", est.pval[[vListIndex]], sep="")
if (out.est.corr) temp5 <- paste(method, ".", est.corr, sep="")
if (!is.null(temp)) {
# NOTE: Names that are not found in a named vector
# get assigned a missing value (NA) to them.
se <- sqrt(diag(temp$cov))[parms]
p <- temp$parms[parms]
if (out.est.beta) outVec[temp1] <- p
if (out.est.se) outVec[temp2] <- se
if (out.est.test) outVec[temp3] <- p/se
if (out.est.pval) outVec[temp4] <- 2*pnorm(abs(p/se), lower.tail=FALSE)
if (out.est.corr) {
# The order is the same as est.corr
corrs <- out.est$corr.parms
cov <- temp$cov
for (i in 1:length(corrs)) {
outVec[temp5[i]] <- cov[corrs[[i]][1], corrs[[i]][2]]
}
}
} else {
if (out.est.beta) outVec[temp1] <- NA
if (out.est.se) outVec[temp2] <- NA
if (out.est.test) outVec[temp3] <- NA
if (out.est.pval) outVec[temp4] <- NA
if (out.est.corr) outVec[temp5] <- NA
}
} # END: for (method in out.est$method)
} # END: if (out.est.flag)
# Effects
if (effectsFlag) {
names <- effects$out.names
nr <- nrow(names)
lnames <- effects$list.names
for (method in effects$method) {
i <- 0
for (type in effects$type) {
i <- i + 1
# Get the list of effects and standard errors
temp <- paste(method, lnames[i], sep="")
leff <- ret[[temp]]
if (is.null(leff)) next
eff <- leff$logEffects
for (j in 1:nr) {
temp <- paste(method, ".eff", type, ".", names[j, ], sep="")
outVec[temp] <- eff[j, ]
}
# Standard errors
eff <- leff$logEffects.se
for (j in 1:nr) {
temp <- paste(method, ".eff", type, ".", names[j, ], ".se", sep="")
outVec[temp] <- eff[j, ]
}
}
}
} # END: if (effectsFlag)
cat(outVec, file=fid, sep="\t")
cat("\n", file=fid)
0
} # END: outputRow
# Function to check the out.est list
check.out.est <- function(out.est) {
corr <- getListName(out.est, "corr.parms")
flag <- !is.null(corr)
new <- corr
if (flag) {
if (class(corr) != "list") {
corr <- unique(corr)
n <- length(corr)
if (n > 1) {
index <- 1
# Convert into a list
new <- list()
for (i in 1:(n-1)) {
for (j in (i+1):n) {
list[[index]] <- c(corr[i], corr[j])
index <- index + 1
}
}
} else {
new <- NULL
}
}
}
out.est$corr.parms <- new
if (!flag) {
out.est <- default.list(out.est, c("parms", "method", "what"),
list(1, "CML", c("beta", "se")),
checkList=list(NA, c("UML", "CML", "EB"),
c("beta", "se", "test", "pvalue")) )
}
out.est
} # END: check.out.est
# Function to check the UML test list
check.tests.uml <- function(u) {
if (is.null(u)) return(NULL)
if (!is.list(u)) {
if (u) {
return(list(omnibus="WALD", main="WALD", inter="WALD"))
} else {
return(NULL)
}
}
tests <- c("WALD", "LR", "SCORE")
omni <- getListName(u, "omnibus")
main <- getListName(u, "main")
inter <- getListName(u, "inter")
temp <- c(omni, main, inter)
temp <- match(temp, tests)
if (any(is.na(temp))) stop("ERROR with tests.UML")
# Get the correct order
if (is.null(omni)) omni <- "WALD"
if (is.null(main)) main <- "WALD"
if (is.null(inter)) inter <- "WALD"
u <- list(omnibus=omni, main=main, inter=inter)
u
} # END: check.tests.uml
# Function to check that the options are consistent
checkOptions <- function(op) {
# Check for the tests.1df option
tests.1df <- getListName(op, "tests.1df")
if (!is.null(tests.1df)) {
op$out.est <- list(parms=tests.1df, what=c("test", "pvalue"))
op$tests.1df <- NULL
}
# Check tests list
tests <- getListName(op, "tests")
if (!is.null(tests)) {
# List with vectors is also of type vector
if (!is.list(tests)) tests <- list(tests)
}
op$tests <- tests
# Make sure at least 1 test is being done
temp <- op$test.omnibus + op$test.main + op$test.inter +
as.numeric(!is.null(getListName(op, "out.est"))) +
as.numeric(!is.null(getListName(op, "tests")))
if (!temp) op$test.omnibus <- 1
op$effects <- checkEffects(getListName(op, "effects"))
op
} # END: checkOptions
# Function to check the effects list
checkEffects <- function(eff) {
if (is.null(eff)) return(NULL)
eff <- default.list(eff,
c("var", "type", "var.levels", "snp.levels", "method", "var.base"),
list("ERROR", 1, 0, 1, c("UML", "CML", "EB"), 0),
error=c(1, 0, 0, 0, 0, 0))
temp <- eff$method %in% c("UML", "CML", "EB")
if (any(temp)) {
eff$method <- eff$method[temp]
} else {
eff$method <- c("UML", "CML", "EB")
}
eff
} # END: checkEffects
# Function to compute the effects
calcEffects <- function(eff, ret) {
for (method in eff$method) {
temp <- ret[[method]]
if (!is.null(temp)) {
for (type in eff$type) {
if (type == 1) {
name <- paste(method, ".JointEffects", sep="")
} else {
name <- paste(method, ".StratEffects", sep="")
}
ret[[name]] <- try(effects.init(temp$parms, temp$cov, eff$var1,
eff$var2, eff$snp.levels, eff$var.levels,
base1=0, base2=eff$var.base, int.var=eff$int.var,
effects=type, sep1=eff$sep), silent=TRUE)
if (class(ret[[name]]) == "try-error") ret[[name]] <- NULL
}
}
}
ret
} # END: calcEffects
# Function to set up effects list
setupEffects <- function(eff, facVars, X.vars, V.vars) {
if (is.null(eff)) return(NULL)
# Call before phenoData0 is removed
eff$var1 <- getMain.vars()[[1]]
temp <- nchar(eff$var1)
if (substring(eff$var1, temp, temp) == "_") {
eff$sep <- ""
} else {
eff$sep <- "_"
}
# Check if the variable is a main effect
temp <- eff$var %in% X.vars
if (!any(temp)) {
temp <- paste(eff$var, " is not a main effect. No effects will be computed",
sep="")
warning(temp)
return(NULL)
}
eff$var <- eff$var[temp]
nvar <- length(eff$var)
if ((nvar == 1) && (eff$var %in% facVars)) {
temp <- levels(factor(phenoData0[, eff$var]))
# For now, first level is the baseline
baseline <- paste(eff$var, "_", temp[1], sep="")
temp <- temp[-1]
# Get the dummy variable names
eff$var2 <- paste(eff$var, "_", temp, sep="")
flag <- 1
} else {
eff$var2 <- eff$var
flag <- 0
}
# Check if var is also an interaction
temp <- eff$var %in% V.vars
if (all(temp)) {
if (!flag) {
# For continuous var
eff$int.var <- paste(eff$var1, eff$var, sep="")
} else {
eff$int.var <- paste(eff$var1, eff$var2, sep="")
}
} else {
eff$int.var <- NULL
}
# Get the names for the output data set
n1 <- length(eff$snp.levels)
n2 <- length(eff$var.levels)
if (flag) n2 <- n2 + 1
names <- matrix(data=" ", nrow=n1, ncol=n2)
if (!flag) {
for (i in 1:n1) {
names[i,] <- paste(eff$var1, eff$snp.levels[i], "_",
eff$var2, "_", eff$var.levels, sep="")
}
} else {
temp <- c(baseline, eff$var2)
for (i in 1:n1) {
names[i,] <- paste(eff$var1, eff$snp.levels[i], "_",
temp, sep="")
}
}
eff$out.names <- names
# Get the list names
temp <- character(length(eff$type))
i <- 1
for (type in eff$type) {
if (type == 1) {
temp[i] <- ".JointEffects"
} else {
temp[i] <- ".StratEffects"
}
i <- i + 1
}
eff$list.names <- temp
eff
} # END: setupEffects
# Check the input lists
snp.list <- check.snp.list(snp.list)
pheno.list <- default.list(pheno.list,
c("response.var"), list("ERROR"), error=c(1))
# Rename main.vars to X.vars, int.vars to V.vars
pheno.list$X.vars <- pheno.list[["main.vars", exact=TRUE]]
pheno.list$V.vars <- pheno.list[["int.vars", exact=TRUE]]
# Determine if fomulas were passed in
temp <- pheno.list[["X.vars", exact=TRUE]]
if (!is.null(temp)) {
main.form <- ("formula" %in% class(temp))
} else {
main.form <- FALSE
}
temp <- pheno.list[["V.vars", exact=TRUE]]
if (!is.null(temp)) {
int.form <- ("formula" %in% class(temp))
} else {
int.form <- FALSE
}
temp <- pheno.list[["strata.var", exact=TRUE]]
if (!is.null(temp)) {
s.form <- ("formula" %in% class(temp))
} else {
s.form <- FALSE
}
formFlag <- main.form + int.form + s.form
temp.list <- op[["temp.list", exact=TRUE]]
temp.list <- check.temp.list(temp.list)
temp <- c("BFGS", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "IRLS")
op <- default.list(op,
c("id", "print", "optimizer", "snpName", "tests.CML",
"tests.EB", "genetic.model", "geno.counts", "subject.counts",
"test.omnibus", "test.main", "test.inter", "tests.UML",
"allele.cc"),
list(1, 0, "BFGS", "SNP_", 1, 1, 0, 0, 0, 0, 0, 0, 1, 1),
checkList=list(NA, 0:1, temp, NA, 0:1, 0:1, 0:3,
0:1, 0:1, 0:1, 0:1, 0:1, NA, 1:2))
op <- checkOptions(op)
op$tests.UML <- check.tests.uml(getListName(op, "tests.UML"))
# Check cc.var
temp <- pheno.list[["cc.var", exact=TRUE]]
if (is.null(temp)) {
op$subject.counts <- 0
if (op$allele.cc == 2) pheno.list$cc.var <- pheno.list$response.var
}
# Check the out.est list
out.est <- op[["out.est", exact=TRUE]]
if (!is.null(out.est)) {
out.est <- check.out.est(out.est)
out.est.flag <- 1
# Set flags for output
out.est.beta <- ("beta" %in% out.est$what)
out.est.se <- ("se" %in% out.est$what)
out.est.test <- ("test" %in% out.est$what)
out.est.pval <- ("pvalue" %in% out.est$what)
out.est.corr <- !is.null(getListName(out.est, "corr.parms"))
} else {
out.est.flag <- 0
}
# The genetic model is taken care of in logistic.SNP
snp.list$genetic.model <- NULL
# Flag for streamed input
stream <- snp.list[["stream", exact=TRUE]]
# All the variables must be given by variable name (not column number)
if ((is.numeric(pheno.list$response.var)) ||
(is.numeric(pheno.list$strata.var)) ||
(is.numeric(pheno.list$factor.vars)) ||
(is.numeric(pheno.list$id.var)) ) {
stop("ERROR: variables must be specified by name, not column number")
}
if ((!main.form) && (is.numeric(pheno.list$X.vars)))
stop("ERROR: variables must be specified by name, not column number")
if ((!int.form) && (is.numeric(pheno.list$V.vars)))
stop("ERROR: variables must be specified by name, not column number")
# Keep only the variables we need
if (!formFlag) {
temp <- c(pheno.list$response.var, pheno.list$strata.var,
pheno.list$X.vars, pheno.list$V.vars, pheno.list$id.var,
pheno.list$cc.var)
pheno.list$keep.vars <- unique(temp)
pheno.list$remove.vars <- NULL
} else {
# Do not remove variables
pheno.list$keep.vars <- NULL
}
pheno.list$remove.miss <- 1
pheno.list$make.dummy <- 0
# Get the data vector of snps
tlist <- list(include.row1=0, include.snps=0, return.type=1, MAF=1,
missing=1, snpNames=1, orderByPheno=1, return.pheno=1)
if (stream) {
tlist$file.index <- 1
temp <- paste("TMP_", temp.list$id, sep="")
tempfile <- getTempfile(temp.list$dir, temp)
tlist$outfile <- tempfile
}
temp <- try(getData.1(snp.list, pheno.list, temp.list, op=tlist),
silent=TRUE)
if (class(temp) == "try-error") {
print(temp)
stop("ERROR loading data")
}
snpData <- temp$data
missing <- temp$missing
snpNames <- temp$snpNames
delimiter <- getDelimiter(snp.list, output=1)
design.X <- NULL
design.V <- NULL
nsnps <- length(snpData)
maf.vec <- temp[["MAF", exact=TRUE]]
maf.flag <- !is.null(maf.vec)
controls <- temp[["controls", exact=TRUE]]
alleles <- temp[["alleles", exact=TRUE]]
allFlag <- !is.null(alleles)
ProbG1 <- temp[["ProbG1", exact=TRUE]]
ProbG1.flag <- !is.null(ProbG1)
# Get the phenotype data
phenoData.list <- temp$phenoData.list
phenoData0 <- phenoData.list$data
if (stream) {
snpFid <- temp$fid
in.miss <- getListName(snp.list, "in.miss")
heter.codes <- getListName(snp.list, "heter.codes")
delete <- temp.list$delete
# Update the snp list
snp.list <- update.snp.list(snp.list, where=1)
snpNames <- getListName(snp.list, "snpNames")
snpFlag <- 1 - is.null(snpNames)
# Get the required objects
temp <- setUp.stream(temp, snp.list, tempfile, delete, controls=controls)
subjFlag <- temp$subjFlag
subj.ids <- temp$subj.ids
total.nsubjects <- temp$total.nsubjects
subj.order2 <- temp$subj.order2
start.stream <- temp$start.stream
stop.stream <- temp$stop.stream
snp <- temp$snp
control.ids <- getListName(temp, "control.ids")
codes <- getInheritanceVec(snp.list$genetic.model,
recode=snp.list$recode)
}
rm(temp, tlist, snp.list, phenoData.list, controls)
temp <- gc(verbose=FALSE)
# Get the response variable
response0 <- as.numeric(phenoData0[, pheno.list$response.var])
if (!any(response0 %in% 0:1)) stop("ERROR: response must be 0-1")
# Get the number of observations
nobs <- length(response0)
# Print out number of observations that will be used
temp <- paste("For the analysis, ", nobs, " observations will be used.", sep="")
print(temp)
print(table(response0))
# Determine all factor variables
facVars <- pheno.list[["factor.vars", exact=TRUE]]
Xflag <- !is.null(pheno.list[["X.vars", exact=TRUE]])
Vflag <- !is.null(pheno.list[["V.vars", exact=TRUE]])
Sflag <- !is.null(pheno.list[["strata.var", exact=TRUE]])
vars <- pheno.list$response.var
if ((!s.form) && (Sflag)) vars <- c(vars, pheno.list$strata.var)
if ((!main.form) && (Xflag)) vars <- c(vars, pheno.list$X.vars)
if ((!int.form) && (Vflag)) vars <- c(vars, pheno.list$V.vars)
vars <- unique(vars)
for (v in vars) {
if ((is.factor(phenoData0[, v])) || (is.character(phenoData0[, v]))) {
facVars <- c(facVars, v)
}
}
facVars <- unique(facVars)
facFlag <- !is.null(facVars)
X.newVars <- NULL
V.newVars <- NULL
# Get the stratafication matrix
temp <- sMatrix.logistic(phenoData0, pheno.list$strata.var, facVars)
design.S0 <- temp$design.S0
op$s.1catVar <- temp$s.1catVar
nStrata <- ncol(design.S0)
# Get the X design matrix
if (Xflag) {
temp <- logistic.dsgnMat(phenoData0, pheno.list$X.vars,
facVars)
design.X0 <- temp$designMatrix
X.newVars <- temp$newVars
} else {
design.X0 <- NULL
}
# Get the V design matrix
if (Vflag) {
temp <- logistic.dsgnMat(phenoData0, pheno.list$V.vars,
facVars)
design.V0 <- temp$designMatrix
V.newVars <- temp$newVars
} else {
design.V0 <- NULL
op$test.inter <- 0
}
# Get the genetic model
genetic.model <- op$genetic.model
# Check for effects
op$effects <- setupEffects(getListName(op, "effects"), facVars,
pheno.list$X.vars, pheno.list$V.vars)
effects <- getListName(op, "effects")
effectsFlag <- !is.null(effects)
# Get a logical vector for obtaining the cc counts
if (op$subject.counts) {
cntrl <- (pheno.list$cc.var == 0)
subj.cnts0 <- table(cntrl)
}
rm(pheno.list, phenoData0)
temp <- gc(verbose=FALSE)
# Set the output directory
if (!is.null(op$out.dir)) {
out.dir <- checkForSep(op$out.dir)
if (nsnps != length(unique(snpNames))) stop("SNP names are not unique")
} else {
out.dir <- NULL
}
# Determine if tests will be computed
omniFlag <- op$test.omnibus
mainFlag <- op$test.main
interFlag <- op$test.inter
test2Flag <- (!is.null(getListName(op, "tests")))
tests.UML <- (!is.null(getListName(op, "tests.UML")))
if (tests.UML) tests.UML.vec <- op$tests.UML
tests.CML <- op$tests.CML
tests.EB <- op$tests.EB
tests <- omniFlag + mainFlag + test2Flag + tests.UML +
tests.CML + tests.EB + interFlag
geno.counts <- op$geno.counts
subj.counts <- op$subject.counts
# Check out.est
if (out.est.flag) {
if ((out.est$parms == 2) && (!Vflag)) out.est$parms <- 1
}
# Get the variable names that will be used
log.vnames <- logistic.vnames(design.X0, design.V0, nStrata,
snpName=op$snpName, out.est=out.est,
genetic.model=genetic.model)
est.p <- log.vnames[["est.p", exact=TRUE]]
est.parms <- log.vnames[["est.parms", exact=TRUE]]
# Get the variable names for the tests
if (mainFlag) main.vars <- getMain.vars()
if (omniFlag) omni.vars <- getOmni.vars()
if (interFlag) inter.vars <- getInter.vars()
if (test2Flag) {
# test2.vars is a list
test2.vars <- getTest2.vars()
op$test2.vars <- NULL
n.tests <- length(test2.vars)
}
rm(facFlag, facVars, X.newVars, V.newVars)
# Open the output file
fileFlag <- 1 - is.null(op$out.file)
if (fileFlag) {
temp <- initOutfile()
fid <- temp$fid
outVec <- temp$outVec
est.se <- temp$est.se
est.test <- temp$est.test
est.pval <- temp$est.pval
est.corr <- temp$est.corr
}
rm(log.vnames)
# Base model
temp <- getListName(op, "base.outfile")
if (!is.null(temp)) {
i <- callGLM(response0, X.main=design.X0, X.int=design.V0, int.vec=NULL,
family=binomial(), prefix="SNP_", retX=TRUE,
retY=TRUE, inc.int.vec=0, int.vec.base=0)
save(i, file=temp)
print(summary(i))
}
op$errorCheck <- 0
ProbG1.vec <- NULL
if (ProbG1.flag) op$imputed <- 1
# Loop over each SNP
i <- 0
while (1) {
if (fileFlag) outVec[] <- NA
i <- i + 1
if (!stream) {
if (i > nsnps) break
snp <- as.numeric(getVecFromStr(snpData[i], delimiter=delimiter))
snp.miss <- missing[i]
snp.name <- snpNames[i]
if (allFlag) {
majMin <- alleles[i]
} else {
majMin <- " "
}
# Get the MAF
if (maf.flag) {
maf <- maf.vec[i]
} else {
maf <- getMAF(snp)
}
# For imputed data
if (ProbG1.flag) ProbG1.vec <- as.numeric(getVecFromStr(ProbG1[i], delimiter=delimiter))
} else {
if (i > 1) {
start.stream <- start.stream + 1
snp <- getNextObs(i, snpFid, snpFlag, snpNames, tempfile, delete,
start.stream, stop.stream, delimiter)
if (is.null(snp)) break
}
snp.name <- snp[1]
if (snpFlag) snpNames <- update.snpNames(snpNames, snp.name)
# Get the correct order and recode
temp <- orderSNP(snp[-1], snp.name, subj.order2, total.nsubjects,
in.miss, heter.codes, subjFlag=subjFlag, subj.ids=subj.ids,
out.genotypes=codes, control.ids=control.ids)
snp <- as.numeric(temp$SNP)
maf <- temp$MAF
snp.miss <- any(is.na(snp))
if (allFlag) {
majMin <- temp$alleles
} else {
majMin <- " "
}
}
# Remove missing values
if (snp.miss) {
temp <- as.logical(!is.na(snp))
if (subj.counts) subj.notMiss <- temp
response <- response0[temp]
snp <- snp[temp]
design.S <- removeOrKeepRows(design.S0, temp, which=1)
if (Xflag) design.X <- removeOrKeepRows(design.X0, temp, which=1)
if (Vflag) design.V <- removeOrKeepRows(design.V0, temp, which=1)
} else {
response <- response0
design.S <- design.S0
design.X <- design.X0
design.V <- design.V0
subj.notMiss <- NULL
}
# Get the genotype counts
genoCounts <- getGenoCounts(snp)
n.genos <- sum(genoCounts > 0)
# Make sure that there is more than 1 genotype
if (n.genos < 2) {
# Do not call snp.logistic
ret <- list()
} else {
# Set options
if (n.genos == 2) {
op$geno.binary <- 1
op$genetic.model <- 0
} else {
op$geno.binary <- 0
op$genetic.model <- genetic.model
}
# Fit the model
ret <- try(snp.main(response, snp, X.main=design.X, X.int=design.V,
X.strata=design.S, ProbG1=ProbG1.vec, op=op), silent=TRUE)
if ("try-error" %in% class(ret)) ret <- list()
}
# Get the list index for test variable names
vListIndex <- getVListIndex(n.genos)
# Add the snp name
ret$snp <- snp.name
# Compute tests
if (tests) {
if (tests.UML) ret <- addToList(ret, "UML")
if (tests.CML) ret <- addToList(ret, "CML")
if (tests.EB) ret <- addToList(ret, "EB")
}
# Effects
if (effectsFlag) ret <- calcEffects(effects, ret)
# Print the return list
if (op$print) print(ret)
# Save the output
if (!is.null(out.dir)) {
temp <- paste(out.dir, "out_", snp.name, ".rda", sep="")
save(ret, file=temp)
}
if (fileFlag) temp <- outputRow()
} # END: while(1)
if (fileFlag) close(fid)
ret
} # END: snp.scan.logistic
# Function to permform a SNP by environment interaction analysis for 1 SNP.
snp.logistic <- function(data, response.var, snp.var, main.vars=NULL,
int.vars=NULL, strata.var=NULL, op=NULL) {
# INPUT:
# data Data frame containing all the data.
# No default
# response.var Name of the binary response variable coded
# as 0 (controls) and 1 (cases)
# No default.
# snp.var Name of the genotype variable coded as 0, 1, 2 (or 0, 1)
# No default.
# main.vars Character vector of the main effects variables or
# a formula.
# The default is NULL
# int.vars Names of all covariates of interest that will
# interact with the SNP variable or a formula.
# The default is NULL
# strata.var Name of the stratification variable(s) or a formula.
# Any character variable or factor will be considered
# categorical. Any numeric variable will be considered
# continuous. An intercept column will be automatically
# included.
# The default is NULL (1 stratum)
# ProbG1.var Variable for Prob(G = 1) or NULL. Not needed if snp.var
# is of length 3.
# The default is NULL.
#####################################################################
# op List with the following names.
# genetic.model 0-2
# 0: trend
# 1: dominant
# 2: recessive
# 3: general
# The default is 0.
# This option has no effect if the snp is binary.
# reltol Stopping tolerance
# The default is 1e-8
# maxiter Maximum number of iterations
# The default is 100
# optimizer One of :"Nelder-Mead", "BFGS", "CG", "L-BFGS-B",
# "SANN", "IRLS"
# The default is "BFGS"
# snpName Name to be used in the SNP variable and interaction
# variables.
# The default is "SNP_".
# debug 0 or 1 to show the IRLS iterations.
# The default is 0
# fit.null 0 or 1 to fit a NULL model which excludes the snp,
# interactions with the snp, and the interacting covariates
# This option takes precedence over zero.vars.
# The default is 0.
# zero.vars Character vector of variable names to fix or a list
# with names "snp.var", "main.vars", int.vars", and
# "strata.var".
# If zero.vars is a list, then the names is the list can
# either be formulas or character vectors.
# The default is 0.
######################################################################
# Check for errors
if (length(response.var) != 1) stop("response.var must be a single variable")
if (!is.data.frame(data)) stop("data must be a data frame")
ProbG1.var <- NULL
n.snp.var <- length(snp.var)
#if (!(n.snp.var %in% c(1, 3))) {
# stop("snp.var must be a single variable or 3 variables for imputed genotypes")
#}
if (n.snp.var != 1) stop("snp.var must be a single variable")
# Check variable names
vlist <- list(response.var=response.var, snp.var=snp.var, main.vars=main.vars,
int.vars=int.vars, strata.var=strata.var, ProbG1.var=ProbG1.var)
vars <- getAllVars(vlist, names=names(vlist))
temp <- !(vars %in% colnames(data))
if (any(temp)) {
print(vars[temp])
stop("The above variables were not found in the input data")
}
# Check if snp.var is in main.vars or int.vars
if (any(snp.var %in% getAllVars(main.vars))) stop("ERROR: main.vars must not contain snp.var")
if (any(snp.var %in% getAllVars(int.vars))) stop("ERROR: int.vars must not contain snp.var")
main.call <- main.vars
int.call <- int.vars
strata.call <- strata.var
op <- default.list(op, c("snpName", "fit.null", "imputed", "genetic.model"),
list("SNP_", 0, 0, 0))
op$imputed <- 0
if ((!is.numeric(snp.var)) && (n.snp.var == 1)) op$snpName <- snp.var
zeroFlag <- 0
zero.vars <- NULL
#if ((n.snp.var > 1) || (!is.null(ProbG1.var))) {
# op$imputed <- 1
#}
if (n.snp.var == 1) {
snp <- as.numeric(unfactor(data[, snp.var]))
snp <- snp[is.finite(snp)]
#if (!all(snp %in% 0:2)) op$imputed <- 1
if (!all(snp %in% 0:2)) stop("snp.var must be coded 0-1-2")
}
imputed <- op$imputed
genetic.model <- op$genetic.model
if (!(genetic.model %in% 0:3)) stop("op$genetic.model must be 0-3")
#if ((imputed) && (genetic.model == 3)) stop("op$genetic.model must be 0-2 for imputed genotypes")
#if ((imputed) && (n.snp.var == 1) && (genetic.model != 0)) {
# stop("op$genetic.model must be 0 for imputed genotypes and length(snp.var) = 1")
#}
#if ((imputed) && (n.snp.var == 1) && (is.null(ProbG1.var))) {
# stop("ProbG1.var must be specified for imputed genotypes and length(snp.var) = 1")
#}
#if ((imputed) && (length(ProbG1.var) > 1)) stop("ProbG1.var must be NULL or of length 1")
main.form <- ("formula" %in% class(main.vars))
int.form <- ("formula" %in% class(int.vars))
s.form <- ("formula" %in% class(strata.var))
ProbG1 <- NULL
# For impute snp with 3 vars, if all 0, then missing
#if (n.snp.var > 1) {
# snp1 <- as.numeric(unfactor(data[, snp.var[1]]))
# snp2 <- as.numeric(unfactor(data[, snp.var[2]]))
# snp3 <- as.numeric(unfactor(data[, snp.var[3]]))
# temp <- (snp1 %in% 0) & (snp2 %in% 0) & (snp3 %in% 0)
# if (any(temp)) data[temp, snp.var] <- NA
#}
# Remove missing values
temp <- getFormulas(vlist)
miss <- c(NA, NaN, Inf, -Inf)
if (length(temp)) data <- applyFormulas(data, temp, remove=miss)
data <- removeMiss.vars(data, vars=vars, miss=miss)
rm(vlist, vars)
temp <- gc(verbose=FALSE)
# Get the response variable
D <- unfactor(data[, response.var])
nobs <- length(D)
# Get Prob(G=1)
#if ((imputed) && (!is.null(ProbG1.var))) ProbG1 <- as.numeric(unfactor(data[, ProbG1.var]))
# Get the snp variable(s)
if (n.snp.var == 1) {
snp <- as.numeric(unfactor(data[, snp.var]))
} else {
# snp1 <- as.numeric(unfactor(data[, snp.var[1]]))
# snp2 <- as.numeric(unfactor(data[, snp.var[2]]))
# snp3 <- as.numeric(unfactor(data[, snp.var[3]]))
# if (genetic.model == 0) {
# snp <- snp2 + 2*snp3
# } else if (genetic.model == 1) {
# snp <- snp2 + snp3
# } else if (genetic.model == 2) {
# snp <- snp3
# }
# Add it ot the data frame
# data[, op$snpName] <- snp
# Save Prob(G = 1) if needed
# if (is.null(ProbG1)) ProbG1 <- snp2
# rm(snp1, snp2, snp3)
# gc()
}
#if (imputed) op$genetic.model <- 0
facVars <- NULL
sflag <- !is.null(strata.var)
# Check for constant strata variable
if ((sflag) && (!s.form) && (length(strata.var) == 1)) {
temp <- unique(makeVector(data[, strata.var]))
if (length(temp) == 1) sflag <- FALSE
}
# Determine if any strata vars are character
if ((sflag) && (!s.form)) {
for (v in strata.var) {
if (is.character(data[, v])) data[, v] <- factor(data[, v])
}
}
# Get the variables that are factors
for (temp in colnames(data)) {
if (is.factor(data[, temp])) facVars <- c(facVars, temp)
}
# Get the V design matrix
design.V0 <- logistic.dsgnMat(data, int.vars, facVars)$designMatrix
int.vars <- colnames(design.V0)
# For fitting a null model
if (op$fit.null) {
# Remove interaction vars, snp
temp <- getAllVars(int.call)
temp <- unique(temp, int.vars)
zero.vars <- temp
int.vars <- NULL
design.V0 <- NULL
zeroFlag <- 1
op$fixParms <- list(parms=snp.var, values=0)
op$zero.vars <- NULL
}
# Get the stratafication matrix
temp <- sMatrix.logistic(data, strata.var, facVars)
design.S0 <- temp$design.S0
op$s.1catVar <- temp$s.1catVar
# Get the X design matrix
design.X0 <- logistic.dsgnMat(data, main.vars, facVars, remove.vars=zero.vars)$designMatrix
# For the zero.vars option
zero.vars <- op[["zero.vars", exact=TRUE]]
if (!is.null(zero.vars)) {
temp <- list(main.vars=design.X0, int.vars=design.V0, strata.var=design.S0)
temp <- apply_zero.vars(zero.vars, temp, snp.var, facVars, data)
op$fixParms <- temp[["fixParms", exact=TRUE]]
temp <- temp$mat.list
design.V0 <- temp[["int.vars", exact=TRUE]]
int.vars <- colnames(design.V0)
design.X0 <- temp[["main.vars", exact=TRUE]]
design.S0 <- temp[["strata.var", exact=TRUE]]
}
main.vars <- colnames(design.X0)
strata.var <- colnames(design.S0)
# Check for non-dummy continuous variables ant print a warning
wflag <- 0
if (!all(design.X0 %in% 0:1)) wflag <- 1
if ((!wflag) && (!all(design.V0 %in% 0:1))) wflag <- 1
if ((!wflag) && (!all(design.S0 %in% 0:1))) wflag <- 1
if (wflag) warning("For stability of the algorithm, continuous variables should be scaled")
# Call the core function
ret <- snp.main(D, snp, X.main=design.X0, X.int=design.V0,
X.strata=design.S0, ProbG1=ProbG1, op=op)
# Add model info
model <- list(data=data, response.var=response.var, snp.var=op$snpName,
main.vars=main.vars, int.vars=int.vars,
strata.var=strata.var, factors=facVars,
snpName=op$snpName, main.call=main.call,
int.call=int.call, strata.call=strata.call)
ret$model.info <- model
ret
} # END: snp.logistic
# Function to permform a SNP by environment interaction analysis for 1 SNP.
snp.main <- function(D, snp, X.main=NULL, X.int=NULL,
X.strata=NULL, ProbG1=NULL, op=NULL) {
# INPUT:
# D Binary response vector coded as 0 (controls) and 1 (cases)
# No default.
# snp Vector of genotypes or 3 column matrix for imputed snp
# No default.
# X.main Design matrix for all covariates of interest, excluding
# the SNP variable. This matrix should NOT include an
# intercept column.
# The default is NULL
# X.int Design matrix for all covariates of interest that will
# interact with the SNP variable.
# This matrix should NOT include an intercept column.
# The default is NULL
# X.strata NULL or a design matrix for the stratification.
# The default is NULL (1 stratum)
# ProbG1 Vector of Prob(G=1) for imputed SNPs
# The default is NULL
#####################################################################
# op List with the following names.
# genetic.model 0-3
# 0: trend
# 1: dominant
# 2: recessive
# 3: general
# The default is 0.
# This option has no effect if the snp is binary.
# reltol Stopping tolerance
# The default is 1e-8
# maxiter Maximum number of iterations
# The default is 100
# optimizer One of :"Nelder-Mead", "BFGS", "CG", "L-BFGS-B",
# "SANN"
# The default is "BFGS"
# snpName Name to be used in the SNP variable and interaction
# variables.
# The default is "SNP_".
# debug 0 or 1 to show the IRLS iterations.
# The default is 0
# errorCheck 0 or 1 to check for errors with the input arguments
# The default is 1
# use.C.code 0 or 1 to call the C code for the BFGS optimizer
# The default is 1
# fit.null 0 or 1 for fitting a null model. All the necessary
# variables should already be removed from the design
# matrices.
# The default is 0
# num.deriv 0 or 1 to use numerical derivatives.
# Numerical derivatives will be used if genetic.model > 0
# The default is 0
# imputed 0 or 1 for imputed SNPs
# s.1catVar 0 or 1 for 1 categorical strata var
# The default is 0
###################################################################
# fixParms List with names "parms" and "value".
# The default is NULL
# parms Vector of parm names to fix
# No default
# values Numeric vector of fixed values
# No default
######################################################################
# RETURN:
# A list with sublists named UML (standard logistic regression),
# CML (conditional ML), and EB (empirical bayes).
# Each sublist contains the parameter estimates and covariance
# matrix. The lists UML and CML also contain the log-likelihood.
# Local function to return the initial estimates
getInit <- function() {
# Define the formula
ff <- "D ~ 1"
if (nx) ff <- paste(ff, " + X.main ", sep="")
if (!fit.null) {
if (!zeroSNP) ff <- paste(ff, "+ fsnp ", sep="")
if (nv) {
if (!gmodel3) {
ff <- paste(ff, "+ fsnp:X.int", sep="")
} else {
ff <- paste(ff, "+ fsnp[,1]:X.int + fsnp[,2]:X.int", sep="")
}
}
}
ff <- as.formula(ff)
# Get variable names for Z
v <- logistic.vnames(X.main, X.int, nStrata, snpName=op$snpName,
genetic.model=genetic.model, fit.null=fit.null, zeroSNP=zeroSNP)
# Call after snp was multiplied by V
fit <- glm(ff, family=binomial(), model=FALSE, x=FALSE, y=FALSE)
# Check the convergence
if (!fit$converged) return(NULL)
UML.parms <- fit$coefficients
loglike <- getLoglike.glm(fit)
cov <- summary(fit)$cov.scaled
cnames <- v$UML.names
#all <- v$all.names
names(UML.parms) <- cnames
parms <- UML.parms
if (zeroSNP) {
# Add the snp back in
pos <- nx + 1
len <- length(parms)
temp <- c(parms[1:pos], 0)
temp2 <- c(cnames[1:pos], v$snp)
if (nv) {
temp <- c(temp, parms[(pos+1):len])
temp2 <- c(temp2, cnames[(pos+1):len])
}
parms <- temp
names(parms) <- temp2
}
# Check for NAs
naFlag <- 0
naXPos <- NULL
naVPos <- NULL
temp <- is.na(UML.parms)
cnames <- cnames[!temp]
colnames(cov) <- cnames
rownames(cov) <- cnames
if (any(temp)) {
naFlag <- 1
# Stop if snp was not estimated
if (temp[nx+2]) stop("ERROR: SNP was not estimated for the UML method")
# Get the NA positions in the design matrices
naPos <- (1:length(UML.parms))[temp]
temp2 <- (naPos <= nx + 1)
if (any(temp2)) naXPos <- naPos[temp2] - 1
temp2 <- !temp2
if (any(temp2)) naVPos <- naPos[temp2] - nx - 2
UML.parms <- UML.parms[!temp]
parms <- parms[!is.na(parms)]
}
# Remove the special character "`" that glm will sometimes
# put in the variable names
parms <- changeStr.names(parms, "`", replace="")
cov <- changeStr.names(cov, "`", replace="")
alpha <- parms[1]
beta <- parms[-1]
# Vector for xi parms. Only intercept will be non zero
xi <- rep(0, times=nStrata)
total <- 2*length(D)
if (op$s.1catVar) {
mini <- NULL
maxi <- NULL
for (i in 1:ncol(X.strata)) {
temp <- X.strata[, i] == 1
freq <- sum(snp[temp])/total
if (freq == 0) {
mini <- c(mini, i)
} else if (freq == 1) {
maxi <- c(maxi, i)
} else {
xi[i] <- log(freq/(1-freq))
}
if (!is.finite(xi[i])) xi[i] <- 0
if (!is.null(mini)) xi[mini] <- min(xi, -5)
if (!is.null(maxi)) xi[maxi] <- max(xi, 5)
}
} else {
freq <- sum(snp)/total
if (freq == 0) {
xi[1] <- -1
} else if (freq == 1) {
xi[1] <- 1
} else {
xi[1] <- log(freq/(1-freq))
}
if (!is.finite(xi[1])) xi[1] <- 0
}
eta <- c(parms, xi)
names(eta) <- c(names(parms), v$strata)
list(eta=eta, alpha=alpha, beta=beta, xi=xi, fit=fit, parms=parms,
loglike=loglike, cov=cov, fitted.values=fit$fitted.values, UML.parms=UML.parms,
naFlag=naFlag, naXPos=naXPos, naVPos=naVPos)
} # END: getInit
# Function to compute Pdg.xs = P(D=d, G=g | X, S)
Pdg.xs <- function(ret, alpha, beta, xi) {
# Get the xi parameters for each observation
if (nStrata == 1) {
temp.xi <- xi
} else {
dim(xi) <- c(nStrata, 1)
temp.xi <- X.strata %*% xi
}
# Make sure that beta is a column vector
dim(beta) <- c(nbeta, 1)
# Get theta for each d, g combination
ret[, d0g0.col] <- 0
ret[, d0g1.col] <- log2 + temp.xi
ret[, d0g2.col] <- 2*temp.xi
ret[, d1g0.col] <- alpha + (Z0 %*% beta)
ret[, d1g1.col] <- alpha + (Z1 %*% beta) + log2 + temp.xi
ret[, d1g2.col] <- alpha + (Z2 %*% beta) + 2*temp.xi
ret <- exp(ret)
# For a binary snp that the user input
if (geno.binary) ret[, g2.col] <- 0
# Compute the sum over (d,g)
sum <- rowSums(ret)
# Divide by sum
ret <- matrixDivideVec(ret, sum)
list(Pdg.matrix=ret, Pdg.rowSums=sum)
} # END: Pdg.xs
# Function to compute P(D=1 | X,S)
Pd1.xs <- function(pmat) {
rowSums(pmat[, d1.col])
} # END: Pd1.xs
# Function to compute E(DG | X,S)
Edg.xs <- function(pmat) {
if (gmodel3) {
return(pmat[, c(d1g1.col, d1g2.col)])
} else if (geno.binary) {
return(pmat[, d1g1.col])
} else {
return(pmat[, d1g1.col] + 2*pmat[, d1g2.col])
}
} # END: Edg.xs
# Function to compute E(G | X,S)
Eg.xs <- function(pmat) {
temp <- pmat
temp[, g2.col] <- 2*temp[, g2.col]
temp <- temp[, c(g1.col, g2.col)]
rowSums(temp)
} # END: Eg.xs
# Function to compute E(DG^2 | X,S) = E(D^2G^2 | X,S)
Edgg.xs <- function(pmat) {
if (geno.binary) {
return(pmat[, d1g1.col])
} else {
return(pmat[, d1g1.col] + 4*pmat[, d1g2.col])
}
} # END: Edgg.xs
# Function to compute E(G^2 | X,S)
Egg.xs <- function(pmat) {
temp <- pmat
temp[, g2.col] <- 4*temp[, g2.col]
temp <- temp[, c(g1.col, g2.col)]
rowSums(temp)
} # END: Egg.xs
# Function to return a logical matrix to make the computation of
# the log-likelihood easier.
getLoglike.mat <- function() {
ret <- matrix(data=FALSE, nrow=n, ncol=nlevels)
if (gmodel3) {
col <- 3*D + 1 + snp
} else {
col <- 3*D + 1 + fsnp
}
for (i in 1:n) ret[i, col[i]] <- TRUE
ret
} # END: getLoglike.mat
# Function to compute the log-likelihood (or testing)
loc.getLoglike <- function(eta) {
alpha <- eta[alpha.row]
beta <- eta[beta.row]
xi <- eta[xi.row]
if (fixFlag) eta <- fixGetEta(eta)
Pdg <- Pdg.xs(Pdg, eta[alpha.row], eta[beta.row], eta[xi.row])
if (!imputed) {
Pdg <- Pdg$Pdg.matrix
ret <- sum(log(Pdg[loglike.mat]))
} else {
# Get the xi parameters for each observation
if (nStrata == 1) {
temp.xi <- xi
} else {
dim(xi) <- c(nStrata, 1)
temp.xi <- X.strata %*% xi
}
# Make sure that beta is a column vector
dim(beta) <- c(nbeta, 1)
vec <- D*(alpha + Z.imp %*% beta) + fsnp*temp.xi + log2*ProbG1
vec <- exp(vec)/Pdg$Pdg.rowSums
ret <- sum(log(vec))
}
ret
} # END: loc.getLoglike
# Function to compute a Z matrix
getZ <- function(gvalue) {
temp.V <- NULL
if (gmodel3) {
temp.snp <- matrix(data=0, nrow=n, ncol=2)
if (gvalue) temp.snp[, gvalue] <- 1
if (nv) {
temp.V <- cbind(matrixMultVec(X.int, temp.snp[, 1]),
matrixMultVec(X.int, temp.snp[, 2]))
}
} else {
temp.snp <- rep(gvalue, times=n)
if (nv) temp.V <- matrixMultVec(X.int, temp.snp)
}
ret <- as.matrix(cbind(X.main, temp.snp, temp.V))
ret
} # END: getZ
# Function to compute Z matrix for imputed snp
getZ.imp <- function(genovec) {
temp.V <- NULL
if (nv) temp.V <- matrixMultVec(X.int, genovec)
ret <- as.matrix(cbind(X.main, genovec, temp.V))
ret
} # END: getZ.imp
# Function to compute W(Y - mu)
getWtYmu <- function(eta=NULL, which=1) {
# eta Vector of parms
# Only needed for which = 1
# which 0 or 1, 1 is for the optimizer
# The default is 1
# Get the matrix of probabilities
if (fixFlag) eta <- fixGetEta(eta)
pmat <- Pdg.xs(Pdg, eta[alpha.row], eta[beta.row], eta[xi.row])$Pdg.matrix
temp1 <- D - Pd1.xs(pmat)
temp2 <- Dsnp - Edg.xs(pmat)
temp3 <- snp - Eg.xs(pmat)
dim(temp1) <- c(1, n)
temp1 <- temp1 %*% X.main
if (gmodel3) {
dim(temp2) <- c(2, n)
} else {
dim(temp2) <- c(1, n)
}
temp2 <- temp2 %*% X.int
dim(temp3) <- c(1, n)
temp3 <- temp3 %*% X.strata
temp <- c(temp1, temp2, temp3)
dim(temp) <- c(nparms, 1)
if (fixFlag) temp <- temp[fixMap]
temp
} # END: getWtYmu
# Function to call the optimizer
callOptim <- function() {
# Set the control list
control <- list(fnscale=-1, maxit=op$maxiter, reltol=op$reltol)
# Call the optimizer
if (op$num.deriv) {
ret <- optim(eta0, loc.getLoglike, method=op$optimizer,
control=control, hessian=TRUE)
} else {
ret <- optim(eta0, loc.getLoglike, gr=getWtYmu, method=op$optimizer,
control=control, hessian=TRUE)
}
# Determine if it converged
conv <- (ret$convergence == 0)
# Get the covariance matrix and score
if (conv) {
conv <- 1
cov <- chol(-ret$hessian)
cov <- chol2inv(cov)
cnames <- names(eta0)
colnames(cov) <- cnames
rownames(cov) <- cnames
} else {
cov <- NA
conv <- 0
}
# Return list
list(parms=ret$par, cov=cov, converged=conv, loglike=ret$value)
} # END: callOptim
# Function to compute empirical bayes estimates
getEB <- function() {
#if (fixFlag) return(NULL)
ids <- c(alpha.row, beta.row)
if (zeroSNP) ids <- ids[-length(ids)]
parm1 <- ret$UML$parms
parm2 <- ret$CML$parms[ids]
cov1 <- ret$UML$cov
vnames <- names(parm1)
dim(parm1) <- NULL
dim(parm2) <- NULL
psi2 <- parm1 - parm2
psi2 <- psi2*psi2
phi <- diag(cov1)
denom <- psi2 + phi
# Compute the new estimates
parms <- parm1*psi2/denom + parm2*phi/denom
temp <- (phi*(phi-psi2))/(denom*denom)
if (length(temp) > 1) {
cmat <- cbind(diag(1-temp), diag(temp))
} else {
cmat <- matrix(c(1-temp, temp), nrow=1)
}
rm(psi2, phi, denom, parm1, parm2)
temp <- gc(verbose=FALSE)
# Set up the Z matrix
temp <- NULL
if (nv) {
temp <- removeOrKeepCols(X.int, 1, which=-1)
if (gmodel3) {
temp <- cbind(matrixMultVec(temp, fsnp[, 1]),
matrixMultVec(temp, fsnp[, 2]))
} else {
temp <- matrixMultVec(temp, fsnp)
}
}
# Define the Z matrix. X has an intercept column as the first column
if (!zeroSNP) {
Z <- cbind(X.main, fsnp, temp)
} else {
Z <- cbind(X.main, temp)
}
# Compute the scores for UML
temp <- D - fitted.values
score1 <- matrixMultVec(Z, temp)
# Get the matrix of probabilities P(D=d, G=g|X,S)
parm2 <- ret$CML$parms
if (fixFlag) parm2 <- fixGetEta(parm2)
pmat <- Pdg.xs(Pdg, parm2[alpha.row], parm2[beta.row], parm2[xi.row])$Pdg.matrix
# Get the scores for the 3 parts of eta (parm2)
temp1 <- D - Pd1.xs(pmat)
temp1 <- matrixMultVec(X.main, temp1)
temp2 <- Dsnp - Edg.xs(pmat)
if (gmodel3) {
temp2 <- cbind(matrixMultVec(X.int, temp2[, 1]),
matrixMultVec(X.int, temp2[, 2]))
} else {
temp2 <- matrixMultVec(X.int, temp2)
}
# If the snp was not in the model as a main effect, then remove the
# first column of temp2.
if (zeroSNP) {
if (ncol(temp2) == 1) {
temp2 <- NULL
} else {
temp2 <- temp2[, -1]
}
}
temp3 <- snp - Eg.xs(pmat)
temp3 <- matrixMultVec(X.strata, temp3)
# Get the scores for CML
if (!fit.null) {
score2 <- cbind(temp1, temp2, temp3)
} else {
score2 <- cbind(temp1, temp3)
}
# Free memory
rm(temp1, temp2, temp3, Z, pmat, parm2)
temp <- gc(verbose=FALSE)
# Compute the covariance of beta.UML and eta
temp <- 0
dim1 <- c(ncol(score1), 1)
dim2 <- c(1, ncol(score2))
for (i in 1:n) {
temp1 <- score1[i, ]
temp2 <- score2[i, ]
dim(temp1) <- dim1
dim(temp2) <- dim2
temp <- temp + (temp1 %*% temp2)
}
cov12 <- cov1 %*% temp %*% ret$CML$cov
# We only need the cov(beta.UML, beta.CML) submatrix
cov12 <- cov12[ids, ids]
rm(score1, score2)
temp <- gc(verbose=FALSE)
# Initialize the covariance matrix
nb <- length(parms)
nb2 <- 2*nb
nbp1 <- nb + 1
cov <- matrix(data=NA, nrow=nb2, ncol=nb2)
cov[1:nb, 1:nb] <- cov1
cov[nbp1:nb2, nbp1:nb2] <- ret$CML$cov[ids, ids]
cov[1:nb, nbp1:nb2] <- cov12
cov[nbp1:nb2, 1:nb] <- t(cov12)
# Return this matrix
UML.CML.cov <- cov12
rownames(UML.CML.cov) <- paste("UML.", vnames, sep="")
colnames(UML.CML.cov) <- paste("CML.", vnames, sep="")
# Obtain the final covariance matrix
cov <- cmat %*% cov %*% t(cmat)
colnames(cov) <- vnames
rownames(cov) <- vnames
list(parms=parms, cov=cov, UML.CML.cov=UML.CML.cov)
} # END: getEB
# Function to call before returning the return list
setReturn <- function(ret) {
class(ret) <- "snp.logistic"
if (!is.null(ret$UML)) {
class(ret$UML) <- "UML"
} else {
return(ret)
}
if (is.null(ret$CML)) return(ret)
if (fixFlag) {
xi.row <- fix_xi.row
if (!xi.row[1]) {
class(ret$CML) <- "CML"
return(ret)
}
}
# Transfrom the strata parms ???
xi <- ret$CML$parms[xi.row]
#exi <- exp(xi)
#ret$CML$strata.parms <- exi/(1 + exi)
ret$CML$strata.parms <- xi
# Remove the strata parms from parms
ret$CML$parms <- ret$CML$parms[-xi.row]
# Get the cov matrix for the strata
xi <- ret$CML$cov[xi.row, xi.row]
dim(xi) <- c(nStrata, nStrata)
# Use the delta method. The derivative is a diagonal matrix
#dexi <- ret$CML$strata.parms/(1 + exi)
#ndexi <- length(dexi)
#if (ndexi > 1) dexi <- diag(dexi)
# Get the covariance matrix
#temp <- dexi %*% xi %*% dexi
#dim(temp) <- c(ndexi, ndexi)
#ret$CML$strata.cov <- temp
ret$CML$strata.cov <- xi
# Set the names
temp <- names(ret$CML$strata.parms)
colnames(ret$CML$strata.cov) <- temp
rownames(ret$CML$strata.cov) <- temp
# Save full cov
ret$CML$cov.full <- ret$CML$cov
# Remove strata parms from cov
temp <- ret$CML$cov[-xi.row, -xi.row]
if (length(temp) == 1) {
dim(temp) <- c(1, 1)
rownames(temp) <- colnames(temp) <- "Intercept"
}
ret$CML$cov <- temp
class(ret$CML) <- "CML"
if (!is.null(ret$EB)) class(ret$EB) <- "EB"
ret
} # END: setReturn
# Function to call glm for a factor snp (for now)
fitGLM <- function() {
snp <- factor(snp)
temp <- callGLM(D, X.main=X.main, X.int=X.int, int.vec=snp,
family=op$family, prefix=op$snpName)
if (temp$converged) {
uml <- list(parms=temp$coefficients)
uml$cov <- summary(temp)$cov.scaled
uml$loglike <- getLoglike.glm(temp)
}
list(UML=uml)
} # END: fitGLM
# Function to get a new eta from fixEta0
fixGetEta <- function(eta) {
ret <- fixEta0
ret[fixMap] <- eta
ret
} # END: fixGetEta
# Wrapper for the C function
CML_EB.R <- function() {
if (op$use.C.code == 0) return(NULL)
if (!is.loaded("CML_EB")) {
warning("CML_EB function is not loaded")
return(NULL)
}
# If initial estimates were passed in, then use them
op_eta0 <- op[["init.parms", exact=TRUE]]
if (!is.null(op_eta0)) {
if (length(op_eta0) != nparms) stop("ERROR: init.parms has the wrong length")
eta0[] <- op_eta0
}
error <- 1
# Make the matrices vectors by row
if (!nx) X.main <- 0
if (!nv) X.int <- 0
X.main <- t(X.main)
dim(X.main) <- NULL
X.int <- t(X.int)
dim(X.int) <- NULL
X.strata <- t(X.strata)
dim(X.strata) <- NULL
# Define the return vector
cml.parms <- double(nparms)
cml.cov <- double(nparms*nparms)
cml.ll <- double(1)
error <- as.integer(1)
nbp1 <- nbeta + 1
eb.parms <- double(nbp1)
eb.cov <- double(nbp1*nbp1)
uml.cov <- ret$UML$cov
uml.parms <- ret$UML$parms
if (zeroSNP) {
# Remove snp and update other variables
temp <- match(op$snpName, names(eta0))
if (!is.na(temp)) {
eta0 <- eta0[-temp]
nparms <- length(eta0)
nbeta <- nbeta - 1
}
}
# Vector for UML-CML cov
uml.cml.cov <- double((nbeta+1)*nparms)
##########################################################################
############## Include PACKAGE="CGEN" when building a package ############
##########################################################################
# Call the C function
temp <- .C("CML_EB", as.double(eta0), as.integer(nparms), as.integer(nbeta),
as.integer(D), as.double(snp), as.integer(n), as.double(X.main), as.integer(nx),
as.double(X.int), as.integer(nv), as.double(X.strata), as.integer(nStrata),
as.integer(genetic.model), as.integer(geno.binary),
as.integer(op$maxiter), as.double(op$reltol), as.integer(op$debug),
as.double(uml.cov), as.double(fitted.values),
as.integer(zeroSNP), as.integer(op$num.deriv), as.integer(imputed), as.double(ProbG1),
as.double(uml.parms), error=error, cml.parms=cml.parms, cml.cov=cml.cov, cml.ll=cml.ll,
eb.parms=eb.parms, eb.cov=eb.cov, uml.cml.cov=uml.cml.cov, PACKAGE="CGEN")
error <- temp$error
if (error) return(list())
# Get the covariance matrix
cml.cov <- matrix(temp$cml.cov, nrow=nparms, byrow=TRUE)
if (any(!is.finite(diag(cml.cov)))) return(list())
cml.parms <- temp$cml.parms
cnames <- names(eta0)
names(cml.parms) <- cnames
colnames(cml.cov) <- cnames
rownames(cml.cov) <- cnames
cnames <- cnames[1:(1+nbeta)]
eb.parms <- temp$eb.parms
names(eb.parms) <- cnames
eb.cov <- matrix(temp$eb.cov, nrow=1+nbeta, byrow=TRUE)
colnames(eb.cov) <- cnames
rownames(eb.cov) <- cnames
# UML-CML matrix
uml.cml.cov <- matrix(temp$uml.cml.cov, byrow=TRUE, ncol=nparms)
rownames(uml.cml.cov) <- paste("UML.", cnames, sep="")
colnames(uml.cml.cov) <- paste("CML.", names(eta0), sep="")
# Return list
list(CML=list(parms=cml.parms, cov=cml.cov, loglike=temp$cml.ll),
EB=list(parms=eb.parms, cov=eb.cov, UML.CML.cov=uml.cml.cov))
} # END: CML_EB.R
# Function to check the initial estimates of CML
check_init <- function() {
new <- eta0
# If initial estimates were passed in, then use them
op_eta0 <- op[["init.parms", exact=TRUE]]
if (!is.null(op_eta0)) {
if (length(op_eta0) != nparms) stop("ERROR: init.parms has the wrong length")
new[] <- op_eta0
}
maxll <- loc.getLoglike(eta0)
maxit <- 100
h <- 0.1
steps <- h*abs(new)
temp <- steps <= 1e-8
steps[temp] <- 0.01
for (i in 1:nparms) {
test <- new
point <- new[i]
step <- steps[i]
flag <- 0
# For 2 directions
for (k in 1:2) {
for (j in 1:maxit) {
point0 <- point + step
test[i] <- point0
ll <- loc.getLoglike(test)
if (ll > maxll) {
maxll <- ll
point <- point0
flag <- 1
} else {
break
}
}
if (flag) {
new[i] <- point
break
} else {
# Try other direction
point <- new[i]
step <- -step
}
}
}
new
} # END: check_init
# Check the options list
op <- default.list(op, c("reltol", "maxiter", "optimizer",
"snpName", "debug", "genetic.model", "errorCheck", "geno.binary",
"use.C.code", "only.UML", "fit.null", "num.deriv", "imputed",
"s.1catVar"),
list(1e-8, 100, "BFGS", "SNP_", 0, 0, 1, 0, 1, 0, 0, 0, 0, 0))
snp.nc <- ncol(snp)
if (is.null(snp.nc)) snp.nc <- 0
if (snp.nc == 3) op$imputed <- 1
op$imputed <- 0 # Changed Mar 11, 2015
fixFlag <- 0
fit.null <- op$fit.null
imputed <- op$imputed
if (imputed) {
if (snp.nc == 3) {
ProbG1 <- snp[, 2]
# Create snp vector
gmodel <- op$genetic.model
if (gmodel == 0) {
snp <- snp[, 2] + 2*snp[, 3]
} else if (gmodel == 1) {
snp <- snp[, 2] + snp[, 3]
} else if (gmodel == 2) {
snp <- snp[, 3]
} else {
stop("Incorrect genetic.model with imputed data")
}
}
op$genetic.model <- 0
if (is.null(ProbG1)) stop("ERROR: ProbG1 is NULL")
} else {
ProbG1 <- 0
}
zeroSNP <- FALSE
if (fit.null) {
X.int <- NULL
op$fixParms <- list(parms=op$snpName, values=0)
zeroSNP <- TRUE
fixFlag <- 1
} else {
fixFlag <- !is.null(op[["fixParms", exact=TRUE]])
if (fixFlag) {
if (op$snpName %in% op$fixParms$parms) zeroSNP <- TRUE
}
}
if (zeroSNP) op$genetic.model <- 0
genetic.model <- op$genetic.model
geno.binary <- 0
if (!(genetic.model %in% c(0, 1, 2, 3))) stop("genetic.model must be 0-3")
# Get the number of genotypes
snp <- unfactor(snp, fun=as.numeric)
usnp <- sort(unique(snp))
n <- length(usnp)
# If the input SNP is binary 0-1, set genetic.model to 0
if (!imputed) {
if (!all(usnp %in% c(0, 1, 2))) stop("snp is not coded correctly")
}
if (n == 1) {
stop("snp only has 1 value")
} else if (n == 2) {
if (genetic.model) warning("genetic.model is set to 0")
genetic.model <- 0
if (all(usnp %in% 0:1)) geno.binary <- 1
}
if (genetic.model %in% 1:2) geno.binary <- 1
# Check D
if (!all(unique(D) %in% c(0, 1))) stop("D is not coded correctly")
if (!is.null(X.strata)) {
if (!is.matrix(X.strata)) stop("X.strata is not a matrix")
}
gmodel3 <- (genetic.model == 3)
n <- length(D)
if (op$optimizer != "BFGS") op$use.C.code <- 0
#if (imputed) dim(ProbG1) <- c(n, 1)
# See if X and V were given
if (is.null(X.main)) {
nx <- 0
} else {
nx <- ncol(X.main)
}
if (is.null(X.int)) {
nv <- 0
} else {
nv <- ncol(X.int)
}
# Get the SNP vector for the specific genetic model
fsnp <- snp
if (genetic.model == 1) {
# Dominant
temp <- snp == 2
fsnp[temp] <- 1
} else if (genetic.model == 2) {
temp <- snp == 1
fsnp[temp] <- 0
temp <- snp == 2
fsnp[temp] <- 1
} else if (genetic.model == 3) {
fsnp <- cbind(as.integer(snp == 1), as.integer(snp == 2))
}
# Check the strata
if (is.null(X.strata)) X.strata <- matrix(data=1, nrow=n, ncol=1)
nStrata <- ncol(X.strata)
# Get the initial estimates
temp <- getInit()
if (is.null(temp)) return(NULL)
# Save the initial estimates and initialize the return list
ret <- list()
ret$UML <- list(parms=temp$UML.parms, cov=temp$cov, loglike=temp$loglike)
if (op$only.UML) setReturn(ret)
nbeta <- length(temp$beta)
eta0 <- temp$eta
alpha.row <- 1
beta.row <- 2:(nbeta+1)
xi.row <- (max(beta.row)+1):length(eta0)
cov <- NULL
nparms <- length(eta0)
# Save the fitted values for empirical bayes
fitted.values <- temp$fitted.values
# If there were NAs in UML, modify the design matrices
if (temp$naFlag) {
if (nx) {
pos <- temp$naXPos
len <- length(pos)
if (len) {
if (len == nx) {
nx <- 0
X.main <- NULL
} else {
X.main <- removeOrKeepCols(X.main, pos, which=-1)
nx <- ncol(X.main)
}
}
}
if (nv) {
pos <- temp$naVPos
len <- length(pos)
if (len) {
if (len == nv) {
nv <- 0
X.int <- NULL
} else {
X.int <- removeOrKeepCols(X.int, pos, which=-1)
nv <- ncol(X.int)
}
}
}
}
# Determine if parameters are to be fixed
if (fixFlag) {
#op$use.C.code <- 0
temp <- op$fixParms
# Define the map for fixed parms
temp2 <- match(temp$parms, names(eta0))
temp2 <- temp2[!is.na(temp2)]
if (!length(temp2)) stop("ERROR with fixParms$parms")
fixEta0 <- eta0
fixEta0[temp2] <- temp$values
fixMap <- (1:nparms)[-temp2]
# Change eta0
eta0 <- eta0[fixMap]
# Get the updated xi.row
temp <- sum(xi.row %in% temp2)
nxi <- length(xi.row) - temp
if (nxi) {
temp <- length(eta0)
fix_xi.row <- (temp-nxi+1):temp
} else {
fix_xi.row <- 0
}
}
if (genetic.model) op$num.deriv <- 1
# Call the C code
clist <- try(CML_EB.R(), silent=TRUE)
if (checkTryError(clist, conv=0)) return(setReturn(ret))
if (!is.null(clist)) {
if (!length(clist)) return(setReturn(ret))
ret$CML <- clist$CML
ret$EB <- clist$EB
return(setReturn(ret))
}
# Initialize
nlevels <- 6
d0g0.col <- 1
d0g1.col <- 2
d0g2.col <- 3
d1g0.col <- 4
d1g1.col <- 5
d1g2.col <- 6
d0.col <- c(d0g0.col, d0g1.col, d0g2.col)
d1.col <- c(d1g0.col, d1g1.col, d1g2.col)
g0.col <- c(d0g0.col, d1g0.col)
g1.col <- c(d0g1.col, d1g1.col)
g2.col <- c(d0g2.col, d1g2.col)
log2 <- log(2)
# Define the Z matrices for efficiency
Z0 <- getZ(0)
Z1 <- getZ(1)
Z2 <- getZ(2)
if (imputed) Z.imp <- getZ.imp(fsnp)
# Initialize the matrix to hold all probabilities P(D=d, G=g| X, S)
Pdg <- matrix(data=0, nrow=n, ncol=nlevels)
# Compute D*snp
if (gmodel3) {
Dsnp <- matrixMultVec(fsnp, D)
} else {
Dsnp <- D*fsnp
}
# Add intercept columns to X and V
X.main <- addIntercept(X.main, nrow=n)
X.int <- addIntercept(X.int, nrow=n)
# Define a logical matrix for the calculation of the log-likelihood
loglike.mat <- getLoglike.mat()
# Update initial estimates if needed
eta0 <- check_init()
temp <- try(callOptim(), silent=TRUE)
if (checkTryError(temp, conv=0)) return(setReturn(ret))
if (!temp$converged) return(setReturn(ret))
temp <- list(parms=temp$parms, cov=temp$cov, loglike=temp$loglike)
ret$CML <- temp
# Empirical Bayes
ret$EB <- getEB()
ret <- setReturn(ret)
ret
} # END: snp.main
# Function to return a vector of variable names
logistic.vnames <- function(X, V, nStrata, snpName="SNP_",
out.est=NULL, genetic.model=0,
fit.null=0, zeroSNP=0) {
# Initialize the return list
ret <- list()
# Check out.est
if (!is.null(out.est)) {
temp <- out.est$parms
if (is.character(temp)) {
ret$est.p <- temp
out.est <- 0
} else {
# Numeric value
out.est <- temp
}
} else {
out.est <- 0
}
save1 <- out.est
ret$int <- "Intercept"
if (!is.null(X)) ret$X <- getVarNames(X, prefix="X")
ret$snp <- getVarNames.snp(prefix=snpName, genetic.model=genetic.model)
if (!is.null(V)) {
temp <- getVarNames.int(V, prefix=snpName, genetic.model=genetic.model, sep=":")
} else {
temp <- NULL
if (out.est == 2) out.est <- 1
}
ret$V <- temp
ret$strata <- paste("Allele_freq.", 1:nStrata, sep="")
ret$all <- c(ret$int, ret$X, ret$snp, ret$V, ret$strata)
if (out.est) {
if (out.est == 3) {
ret$est.p <- c(ret$int, ret$X, ret$snp, ret$V)
} else {
ret$est.p <- ret$snp
if (out.est == 2) {
ret$est.p <- c(ret$est.p, ret$V)
}
}
if (genetic.model == 3) {
# Create another vector of parm names for the case when there is
# less than 3 genotypes
snp1 <- paste(snpName, 1, sep="")
if (!is.null(V)) {
tempv <- getVarNames.int(V, prefix=snpName, genetic.model=0, sep=":")
tempv1 <- getVarNames.int(V, prefix=snp1, genetic.model=0, sep=":")
} else {
tempv <- NULL
tempv1 <- NULL
}
if (save1 == 3) {
temp <- c(ret$int, ret$X, snpName, tempv)
temp1 <- c(ret$int, ret$X, snp1, tempv1)
} else if (save1 == 2) {
temp <- c(snpName, tempv)
temp1 <- c(snp1, tempv1)
} else if (save1 == 1) {
temp <- snpName
temp1 <- snp1
}
ret$est.p <- list(ret$est.p, temp)
ret$est.parms <- list(ret$est.p[[1]], temp1)
} else {
ret$est.parms <- ret$est.p
}
} else {
ret$est.parms <- ret$est.p
}
if (fit.null) {
ret$UML.names <- c(ret$int, ret$X)
ret$all.names <- c(ret$int, ret$X, ret$snp, ret$strata)
} else {
UML.names <- c(ret$int, ret$X)
if (!zeroSNP) UML.names <- c(UML.names, ret$snp)
UML.names <- c(UML.names, ret$V)
ret$UML.names <- UML.names
ret$all.names <- ret$all
}
ret
} # END: logistic.vnames
# Function to create a design matrix
logistic.dsgnMat <- function(data, vars, facVars, removeInt=1,
norm.names=1, remove.vars=NULL) {
# data Data frame
# vars Character vector of variable names or a formula
# facVars Character vector of factor names
# remove.vars Character vector of variable names to remove.
# They are removed after the variable names are normalized.
# The default is NULL.
rm.vars <- function(mat, rmvars) {
mnames <- colnames(mat)
temp <- !(mnames %in% rmvars)
mnames <- mnames[temp]
len <- length(mnames)
if (len == 0) {
mat <- NULL
} else if ((len == 1) && (mnames == "Intercept")) {
mat <- NULL
} else {
mat <- removeOrKeepCols(mat, mnames, which=1)
}
mat
} # END: rm.vars
if (is.null(vars)) return(list(designMatrix=NULL, newVars=NULL))
# Determine if vars is a formula
if ("formula" %in% class(vars)) {
# Get the design matrix
design <- model.matrix(vars, data=data)
# Remove the intercept, if needed
newVars <- colnames(design)
if (removeInt) {
if (newVars[1] == "(Intercept)") {
design <- removeOrKeepCols(design, 1, which=-1)
newVars <- newVars[-1]
}
}
if (norm.names) colnames(design) <- normVarNames(colnames(design))
if (!is.null(remove.vars)) design <- rm.vars(design, remove.vars)
return(list(designMatrix=design, newVars=newVars))
}
design <- removeOrKeepCols(data, vars, which=1)
newVars <- NULL
if (!is.null(facVars)) {
temp <- vars %in% facVars
if (any(temp)) {
temp <- vars[temp]
temp <- createDummy(design, vars=temp)
design <- temp$data
newVars <- temp$newVars
}
}
design <- as.matrix(design)
if (norm.names) colnames(design) <- normVarNames(colnames(design))
# Check for constant variables
design <- checkForConstantVar(design, msg=1)$data
if (!removeInt) {
# Add intercept
cnames <- colnames(design)
design <- cbind(1, design)
colnames(design) <- c("Intercept", cnames)
}
if (!is.null(remove.vars)) design <- rm.vars(design, remove.vars)
# Make sure matrix is numeric
d <- dim(design)
cnames <- colnames(design)
design <- as.numeric(design)
dim(design) <- d
colnames(design) <- cnames
list(designMatrix=design, newVars=newVars)
} # END: logistic.dsgnMat
# Function to apply the zero.vars option to the design matrices
apply_zero.vars <- function(zero.vars, mat.list, snp.var, facVars, data) {
mat.names <- names(mat.list)
fixParms <- NULL
if (is.character(zero.vars)) {
for (i in 1:length(mat.names)) {
mat <- mat.list[[i]]
if (is.null(mat)) next
cnames <- colnames(mat)
temp <- cnames %in% zero.vars
if (any(temp)) {
cnames <- cnames[!temp]
if (!length(cnames)) {
mat <- NULL
} else {
mat <- removeOrKeepCols(mat, cnames, which=1)
}
mat.list[[mat.names[i]]] <- mat
}
}
if (snp.var %in% zero.vars) {
fixParms <- list(parms=snp.var, values=0)
}
} else {
# zero.vars is a list
znames <- names(zero.vars)
temp <- zero.vars[["snp.var", exact=TRUE]]
if (!is.null(temp)) {
if (temp == snp.var) {
fixParms <- list(parms=snp.var, values=0)
}
}
for (z in znames) {
vars <- zero.vars[[z]]
if (is.null(vars)) next
mat <- mat.list[[z, exact=TRUE]]
if (is.null(mat)) next
cnames <- colnames(mat)
pnames <- colnames(logistic.dsgnMat(data, vars, facVars)$designMatrix)
temp <- cnames %in% pnames
if (any(temp)) {
cnames <- cnames[!temp]
if (!length(cnames)) {
mat <- NULL
} else {
mat <- removeOrKeepCols(mat, cnames, which=1)
}
mat.list[[z]] <- mat
}
}
}
list(mat.list=mat.list, fixParms=fixParms)
} # END: apply_zero.vars
# Function to output the needed UML and CML estimates
UML_CML_GxE_parms <- function(main.variables, interaction.variables, out.file, snps, data, response.var,
main.vars=NULL, int.vars=NULL, strata.var=NULL, ProbG1.var=NULL, op=NULL) {
# Input arguments:
# out.file Output file to save the necessary UML and CML estimates, or NULL
# snps Character vector of snps to be analyzed
# data Data frame containing the snps to be analyzed, outcome,
# covariates and other variables
# response.var Same definition as is snp.logistic
# main.vars Same definition as is snp.logistic
# int.var Character vector of interaction variable
# strata.var Same definition as is snp.logistic
# op Same definition as is snp.logistic
# Output variable names for the UML main effects:
# UML.G.BETA UML.E.BETA UML.GE.BETA
# Output variable names for the upper triangle of the UML covariance matrix:
# UML.G.G.COV, UML.G.E.COV UML.G.GE.COV
# UML.E.E.COV UML.E.GE.COV
# UML.GE.GE.COV
# Output variable names for the UML-CML covariance matrix:
# UML.G.CML.G.COV, UML.G.CML.E.COV UML.G.CML.GE.COV
# UML.E.CML.G.COV, UML.E.CML.E.COV UML.E.CML.GE.COV
# UML.GE.CML.G.COV, UML.GE.CML.E.COV UML.GE.CML.GE.COV
n.main.vars <- length(main.variables)
n.int.vars <- length(interaction.variables)
UML.ME.out <- outNames.me("UML", n.main.vars, n.int.vars)
CML.ME.out <- outNames.me("CML", n.main.vars, n.int.vars)
UML.COV.out <- outNames.cov("UML", n.main.vars, n.int.vars)
CML.COV.out <- outNames.cov("CML", n.main.vars, n.int.vars)
UML.CML.COV.out <- outNames.cov2(n.main.vars, n.int.vars)
fileFlag <- !is.null(out.file)
print <- op$print
# Open the output file
if (fileFlag) fid <- file(out.file, "w")
# Write out the column names
outvars <- c("SNP", UML.ME.out, UML.COV.out, CML.ME.out, CML.COV.out, UML.CML.COV.out)
if (fileFlag) writeOut(fid, outvars)
# Vector to store output values
outvec <- character(length(outvars))
names(outvec) <- outvars
outvec0 <- outvec
# Loop over each snp
for (snp in snps) {
outvec0[] <- ""
outvec <- outvec0
outvec["SNP"] <- snp
fit <- try(snp.logistic(data, response.var, snp, main.vars=main.vars,
int.vars=int.vars, strata.var=strata.var, op=op), silent=TRUE)
if ("try-error" %in% class(fit)) {
if (print) print(fit)
if (fileFlag) writeOut(fid, outvec)
next
}
if (print) print(summary(fit))
# Extract UML and CML estimates
temp <- try(extractEst(fit, "UML", outvec, snp, main.variables, interaction.variables,
UML.ME.out, UML.COV.out), silent=TRUE)
if (!("try-error" %in% class(temp))) outvec <- temp
temp <- try(extractEst(fit, "CML", outvec, snp, main.variables, interaction.variables,
CML.ME.out, CML.COV.out), silent=TRUE)
if (!("try-error" %in% class(temp))) outvec <- temp
temp <- try(extract_UML.CML(fit, outvec, snp, main.variables, interaction.variables,
UML.CML.COV.out), silent=TRUE)
if (!("try-error" %in% class(temp))) outvec <- temp
if (fileFlag) writeOut(fid, outvec)
}
if (fileFlag) close(fid)
outvec
} # END: UML_CML_GxE_parms
# Function to write output to an open (tab-delimited) file
writeOut <- function(FID, vec) {
# FID File id
# vec Vector to output
str <- paste(vec, collapse="\t", sep="")
write(str, file=FID, ncolumns=1)
NULL
} # END: writeOut
# Function to get the parameter names
outNames.me <- function(which, n.main.vars, n.int.vars) {
ret <- "G"
if (n.main.vars < 2) {
ret <- c(ret, "E")
} else {
jj <- 1:n.main.vars
ret <- c(ret, paste("E", jj, sep=""))
}
if (n.int.vars < 2) {
ret <- c(ret, "GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("GE", jj, sep=""))
}
ret <- paste(which, ".", ret, ".BETA", sep="")
ret
} # END: outNames.me
# Function to get output names
outNames.cov <- function(which, n.main.vars, n.int.vars) {
# G row
ret <- "G.G"
if (n.main.vars < 2) {
ret <- c(ret, "G.E")
} else {
jj <- 1:n.main.vars
ret <- c(ret, paste("G.E", jj, sep=""))
}
if (n.int.vars < 2) {
ret <- c(ret, "G.GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("G.GE", jj, sep=""))
}
# E rows
if (n.main.vars < 2) {
ret <- c(ret, "E.E")
if (n.int.vars < 2) {
ret <- c(ret, "E.GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("E.GE", jj, sep=""))
}
} else {
jj0 <- 1:n.main.vars
jj2 <- 1:n.int.vars
for (i in 1:n.main.vars) {
jj <- i:n.main.vars
ret <- c(ret, paste("E", i, ".E", jj, sep=""))
if (n.int.vars < 2) {
ret <- c(ret, paste("E", i, ".GE", sep=""))
} else {
ret <- c(ret, paste("E", i, ".GE", jj2, sep=""))
}
}
}
# GE rows
if (n.int.vars < 2) {
ret <- c(ret, "GE.GE")
} else {
jj0 <- 1:n.int.vars
for (i in 1:n.int.vars) {
jj <- i:n.int.vars
ret <- c(ret, paste("GE", i, ".GE", jj, sep=""))
}
}
ret <- paste(which, ".", ret, ".COV", sep="")
ret
} # END: outNames.cov
outNames.cov2 <- function(n.main.vars, n.int.vars) {
if ((n.main.vars < 2) && (n.int.vars < 2)) {
ret <- c("UML.G.CML.G.COV", "UML.G.CML.E.COV", "UML.G.CML.GE.COV",
"UML.E.CML.G.COV", "UML.E.CML.E.COV", "UML.E.CML.GE.COV",
"UML.GE.CML.G.COV", "UML.GE.CML.E.COV", "UML.GE.CML.GE.COV")
return(ret)
}
# G row
ret <- "UML.G.CML.G"
if (n.main.vars < 2) {
ret <- c(ret, "UML.G.CML.E")
} else {
jj <- 1:n.main.vars
ret <- c(ret, paste("UML.G.CML.E", jj, sep=""))
}
if (n.int.vars < 2) {
ret <- c(ret, "UML.G.CML.GE")
} else {
jj <- 1:n.int.vars
ret <- c(ret, paste("UML.G.CML.GE", jj, sep=""))
}
# E rows
if (n.main.vars < 2) {
ret <- c(ret, "UML.E.CML.G")
ret <- c(ret, "UML.E.CML.E")
if (n.int.vars < 2) {
ret <- c(ret, "UML.E.CML.GE")
} else {
ret <- c(ret, paste("UML.E.CML.GE", 1:n.int.vars, sep=""))
}
} else {
jj <- 1:n.main.vars
jj2 <- 1:n.int.vars
for (i in 1:n.main.vars) {
ret <- c(ret, paste("UML.E", i, ".CML.G", sep=""))
ret <- c(ret, paste("UML.E", i, ".CML.E", jj, sep=""))
if (n.int.vars < 2) {
ret <- c(ret, paste("UML.E", i, ".CML.GE", sep=""))
} else {
ret <- c(ret, paste("UML.E", i, ".CML.GE", jj2, sep=""))
}
}
}
# GE rows
if (n.int.vars < 2) {
ret <- c(ret, "UML.GE.CML.G")
if (n.main.vars < 2) {
ret <- c(ret, "UML.GE.CML.E")
} else {
ret <- c(ret, paste("UML.GE.CML.E", 1:n.main.vars, sep=""))
}
ret <- c(ret, "UML.GE.CML.GE")
} else {
jj0 <- 1:n.main.vars
jj2 <- 1:n.int.vars
for (i in 1:n.int.vars) {
ret <- c(ret, paste("UML.GE", i, ".CML.G", sep=""))
if (n.main.vars < 2) {
ret <- c(ret, paste("UML.GE", i, ".CML.E", sep=""))
} else {
ret <- c(ret, paste("UML.GE", i, ".CML.E", jj0, sep=""))
}
ret <- c(ret, paste("UML.GE", i, ".CML.GE", jj2, sep=""))
}
}
ret <- paste(ret, ".COV", sep="")
ret
} # END: outNames.cov2
# Function to get the UML-CML cov estimates
extract_UML.CML <- function(fit, outvec, G.name, E.name, GE.name, ids) {
x <- fit[["EB", exact=TRUE]]
if (is.null(x)) return(outvec)
cov <- x[["UML.CML.cov", exact=TRUE]]
if (is.null(cov)) return(outvec)
# Get the names of the interaction parms
GE.name <- paste(G.name, ":", GE.name, sep="")
G.row <- paste("UML.", G.name, sep="")
G.col <- paste("CML.", G.name, sep="")
E.row <- paste("UML.", E.name, sep="")
E.col <- paste("CML.", E.name, sep="")
GE.row <- paste("UML.", GE.name, sep="")
GE.col <- paste("CML.", GE.name, sep="")
# Vector of expected parameter names
row.names <- paste("UML.", c(G.name, E.name, GE.name), sep="")
col.names <- paste("CML.", c(G.name, E.name, GE.name), sep="")
n.names <- length(col.names)
n.int.vars <- length(E.name)
cnames <- colnames(cov)
temp <- col.names %in% cnames
if (!all(temp)) {
cov2 <- matrix(data=NA, nrow=n.names, ncol=n.names)
colnames(cov2) <- col.names
rownames(cov2) <- row.names
col2.names <- col.names[temp]
row2.names <- row.names[temp]
cov2[row2.names, col2.names] <- cov[row2.names, col2.names]
} else {
cov2 <- cov
}
if (n.int.vars < 2) {
outvec[ids] <- c(cov2[G.row, G.col], cov2[G.row, E.col], cov2[G.row, GE.col],
cov2[E.row, G.col], cov2[E.row, E.col], cov2[E.row, GE.col],
cov2[GE.row, G.col], cov2[GE.row, E.col], cov2[GE.row, GE.col])
} else {
vec <- cov2[G.row, G.col]
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.row, E.col[i]])
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.row, GE.col[i]])
for (i in 1:n.int.vars) {
row <- E.row[i]
vec <- c(vec, cov2[row, G.col])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, E.col[j]])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, GE.col[j]])
}
for (i in 1:n.int.vars) {
row <- GE.row[i]
vec <- c(vec, cov2[row, G.col])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, E.col[j]])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, GE.col[j]])
}
outvec[ids] <- vec
}
outvec
} # END: extract_UML.CML
# Function to extract estimates. The updated output vector will be returned
extractEst <- function(fit, which, outvec, G.name, E.name, GE.name, ids.main, ids.cov) {
# fit Return object from snp.logistic
# which "UML" or "CML"
# outvec Output vector
# G.name Name of the SNP
# E.name Name of the environmental variable
if (which == "EB") return(extract_UML.CML(fit, outvec, G.name, E.name))
x <- fit[[which, exact=TRUE]]
if (is.null(x)) return(outvec)
n.int.vars <- length(E.name)
# Get the names of the interaction parms
GE.name <- paste(G.name, ":", GE.name, sep="")
# Vector of expected parameter names
vec.names <- c(G.name, E.name, GE.name)
n.vec.names <- length(vec.names)
parms <- x[["parms", exact=TRUE]]
if (!is.null(parms)) {
cnames <- names(parms)
temp <- vec.names %in% cnames
if (!all(temp)) {
parms2 <- rep(NA, n.vec.names)
names(parms2) <- vec.names
vec2.names <- vec.names[temp]
parms2[vec2.names] <- parms[vec2.names]
} else {
parms2 <- parms
vec2.names <- vec.names
}
outvec[ids.main] <- c(parms[vec2.names])
}
cov <- x[["cov", exact=TRUE]]
if (!is.null(cov)) {
cnames <- colnames(cov)
temp <- vec.names %in% cnames
if (!all(temp)) {
cov2 <- matrix(data=NA, nrow=n.vec.names, ncol=n.vec.names)
colnames(cov2) <- vec.names
rownames(cov2) <- vec.names
vec2.names <- vec.names[temp]
cov2[vec2.names, vec2.names] <- cov[vec2.names, vec2.names]
} else {
cov2 <- cov
}
if (n.int.vars < 2) {
outvec[ids.cov] <- c(cov2[G.name, G.name], cov2[G.name, E.name], cov2[G.name, GE.name],
cov2[E.name, E.name], cov2[E.name, GE.name], cov2[GE.name, GE.name])
} else {
vec <- cov2[G.name, G.name]
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.name, E.name[i]])
for (i in 1:n.int.vars) vec <- c(vec, cov2[G.name, GE.name[i]])
for (i in 1:n.int.vars) {
row <- E.name[i]
for (j in i:n.int.vars) vec <- c(vec, cov2[row, E.name[j]])
for (j in 1:n.int.vars) vec <- c(vec, cov2[row, GE.name[j]])
}
for (i in 1:n.int.vars) {
row <- GE.name[i]
for (j in i:n.int.vars) vec <- c(vec, cov2[row, GE.name[j]])
}
outvec[ids.cov] <- vec
}
}
outvec
} # END: extractEst
# Function for running a scan
scan.UML_CML <- function(snp.list, pheno.list, op=NULL) {
# Check the input lists
snp.list <- check.snp.list(snp.list)
pheno.list <- default.list(pheno.list,
c("response.var"), list("ERROR"), error=c(1))
response.var <- pheno.list$response.var
main.vars <- pheno.list[["main.vars", exact=TRUE]]
int.vars <- pheno.list[["int.vars", exact=TRUE]]
strata.var <- pheno.list[["strata.vars", exact=TRUE]]
cc.var <- pheno.list[["cc.var", exact=TRUE]]
# Determine if fomulas were passed in
temp <- main.vars
if (!is.null(temp)) {
main.form <- ("formula" %in% class(temp))
} else {
main.form <- FALSE
}
temp <- int.vars
if (!is.null(temp)) {
int.form <- ("formula" %in% class(temp))
} else {
int.form <- FALSE
}
temp <- strata.var
if (!is.null(temp)) {
s.form <- ("formula" %in% class(temp))
} else {
s.form <- FALSE
}
formFlag <- main.form + int.form + s.form
temp.list <- op[["temp.list", exact=TRUE]]
temp.list <- check.temp.list(temp.list)
temp <- c("BFGS", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "IRLS")
op <- default.list(op,
c("id", "print", "optimizer", "snpName","genetic.model",
"allele.cc", "out.file", "save.varnames"),
list(1, 0, "BFGS", "SNP_", 0, 1, "scan.UML_CML.txt", 1),
error=c(0, 0, 0, 0, 0, 0, 0, 0),
checkList=list(NA, 0:1, temp, NA, 0:2, 1:2, NA, 0:1))
print <- op$print
SNP <- op$snpName
# Check cc.var
if (is.null(cc.var)) {
if (op$allele.cc == 2) pheno.list$cc.var <- pheno.list$response.var
}
# The genetic model is taken care of in logistic.SNP
snp.list$genetic.model <- NULL
# All the variables must be given by variable name (not column number)
if ((is.numeric(pheno.list$response.var)) ||
(is.numeric(pheno.list$strata.var)) ||
(is.numeric(pheno.list$factor.vars)) ||
(is.numeric(pheno.list$id.var)) ) {
stop("ERROR: variables must be specified by name, not column number")
}
if ((!main.form) && (is.numeric(pheno.list$X.vars)))
stop("ERROR: variables must be specified by name, not column number")
if ((!int.form) && (is.numeric(pheno.list$V.vars)))
stop("ERROR: variables must be specified by name, not column number")
# Keep only the variables we need
if (!formFlag) {
temp <- c(response.var, strata.var, main.vars, int.vars,
pheno.list$id.var, cc.var)
pheno.list$keep.vars <- unique(temp)
pheno.list$remove.vars <- NULL
} else {
# Do not remove variables
pheno.list$keep.vars <- NULL
}
pheno.list$remove.miss <- 1
pheno.list$make.dummy <- 0
# Get the data vector of snps
tlist <- list(include.row1=0, include.snps=0, return.type=1, MAF=1,
missing=1, snpNames=1, orderByPheno=1, return.pheno=1, ProbG1=1)
temp <- try(getData.1(snp.list, pheno.list, temp.list, op=tlist),
silent=TRUE)
if (class(temp) == "try-error") {
print(temp)
stop("ERROR loading data")
}
snpData <- temp$data
snpNames <- temp$snpNames
delimiter <- getDelimiter(snp.list, output=1)
nsnps <- length(snpData)
maf.vec <- temp[["MAF", exact=TRUE]]
maf.flag <- !is.null(maf.vec)
alleles <- temp[["alleles", exact=TRUE]]
allFlag <- !is.null(alleles)
ProbG1 <- temp[["ProbG1", exact=TRUE]]
ProbG1.flag <- !is.null(ProbG1)
ProbG1.var <- NULL
# Get the phenotype data
phenoData.list <- temp$phenoData.list
phenoData0 <- phenoData.list$data
phenoData0 <- as.data.frame(phenoData0, stringsAsFactors=FALSE)
phenoData0[, SNP] <- NA
if (ProbG1.flag) {
ProbG1.var <- paste(SNP, "ProbG1", sep="")
phenoData0[, ProbG1.var] <- NA
}
rm(temp, tlist, snp.list, phenoData.list)
temp <- gc(verbose=FALSE)
# Get the response variable
response0 <- as.numeric(phenoData0[, response.var])
if (!any(response0 %in% 0:1)) stop("ERROR: response must be 0-1")
# Get the number of observations
nobs <- length(response0)
# Print out number of observations that will be used
temp <- paste("For the analysis, ", nobs, " observations will be used.\n", sep="")
cat(temp)
cat(paste("Number of cases = ", sum(response0 %in% 1), "\n", sep=""))
cat(paste("Number of controls = ", sum(response0 %in% 0), "\n", sep=""))
cat(paste("Output will be written to: ", op$out.file, "\n\n", sep=""))
rm(pheno.list)
temp <- gc(verbose=FALSE)
# Run a base model with simulated SNP
phenoData0[, SNP] <- rbinom(nobs, 1, 0.5)
out.file <- op$out.file
out.base <- paste(out.file, "_info", sep="")
temp <- try(snp.logistic(phenoData0, response.var, SNP, main.vars=main.vars,
int.vars=int.vars, strata.var=NULL, op=op), silent=TRUE)
if ("try-error" %in% class(temp)) {
print(temp)
stop()
}
interaction.variables <- temp$model.info$int.vars
main.variables <- op[["E.parm.names", exact=TRUE]]
if (is.null(main.variables)) main.variables <- interaction.variables
if (op$save.varnames) {
sink(out.base)
cat("Exposure variables:\n")
print(main.variables)
cat("Interaction variables:\n")
print(interaction.variables)
print(summary(temp))
print(table(response0))
for (v in unique(c(main.variables, interaction.variables))) {
cat(paste(v, "\n", sep=""))
tab <- table(phenoData0[, v], response0, exclude=NULL)
if (length(tab) < 31) {
tab <- addmargins(tab)
print(tab)
}
}
sink()
}
if (print) {
cat("Exposure variables:\n")
print(main.variables)
cat("Interaction variables:\n")
print(interaction.variables)
}
rm(response0)
gc()
# Make sure all interaction variables are also main effects
temp <- temp$UML$parms
temp <- names(temp)
if (!(all(main.variables %in% temp))) {
stop("All interaction variables must also be main effects")
}
# Flag for first line of output
firstFlag <- 0
# Open the output file
fid <- file(op$out.file, "w")
# Loop over each SNP
i <- 0
while (1) {
errorFlag <- 0
i <- i + 1
if (i > nsnps) break
snp <- as.numeric(getVecFromStr(snpData[i], delimiter=delimiter))
snp.name <- snpNames[i]
phenoData0[, SNP] <- snp
if (allFlag) {
majMin <- alleles[i]
} else {
majMin <- " "
}
# Get the MAF
if (maf.flag) {
maf <- maf.vec[i]
} else {
maf <- getMAF(snp)
}
# For imputed data
if (ProbG1.flag) {
temp <- as.numeric(getVecFromStr(ProbG1[i], delimiter=delimiter))
phenoData0[, ProbG1.var] <- temp
}
# Fit the model
temp <- try(UML_CML_GxE_parms(main.variables, interaction.variables, NULL, SNP, phenoData0, response.var,
main.vars=main.vars, int.vars=int.vars, strata.var=strata.var, ProbG1.var=ProbG1.var,
op=op), silent=TRUE)
errorFlag <- ("try-error" %in% class(temp))
if (!errorFlag) {
if (firstFlag) {
writeOut(fid, c(snp.name, majMin, maf, temp[-1]))
} else {
writeOut(fid, c("SNP", "Alleles", "MAF", names(temp[-1])))
writeOut(fid, c(snp.name, majMin, maf, temp[-1]))
firstFlag <- 1
errorVec <- temp[-1]
errorVec[] <- "NA"
}
} else {
if (firstFlag) writeOut(fid, c(snp.name, majMin, maf, errorVec))
}
} # END: while(1)
close(fid)
op$out.file
} # END: scan.UML_CML
# Function for meta-anlaysis
meta.fixed <- function(study.list) {
# study.list List of named lists with names beta and cov
nstudy <- length(study.list)
sigma.inv.list <- list()
sum.sigma.inv <- 0
sum.sigma.inv.beta <- 0
# Compute sigma inverse for each study
for (i in 1:nstudy) {
temp <- study.list[[i]]
beta <- temp$beta
sigma <- temp$cov
sigma.inv <- solve(sigma)
sum.sigma.inv <- sum.sigma.inv + sigma.inv
dim(beta) <- c(length(beta), 1)
sum.sigma.inv.beta <- sum.sigma.inv.beta + sum.sigma.inv %*% beta
sigma.inv.list[[i]] <- sigma.inv
}
meta.cov <- solve(sum.sigma.inv)
meta.beta <- meta.cov %*% sum.sigma.inv.beta
list(meta.beta=meta.beta, meta.cov=meta.cov, meta.cov.inv=sum.sigma.inv,
inv.list=sigma.inv.list)
} # END: meta.fixed
# Function to compute A matrix for EB meta
get.EB.A_matrix <- function(parm1, cov1, parm2) {
psi <- parm1 - parm2
dim(psi) <- NULL
psi2 <- psi*psi
phi <- diag(cov1)
temp <- psi2/(psi2 + phi)
A <- diag(temp)
A
} # END: get.EB.A_matrix
meta.EB <- function(beta1, beta2, sumInv1.inv, sumInv2.inv, inv1.list, inv2.list,
cov12.list) {
# Covariance matrix for beta1 is sumInv1.inv
A <- get.EB.A_matrix(beta1, sumInv1.inv, beta2)
n <- length(beta1)
I <- matrix(0, nrow=n, ncol=n)
diag(I) <- 1
beta.EB <- (A %*% beta1) + (I - A) %*% beta2
sum <- 0
n.study <- length(inv1.list)
for (i in 1:n.study) {
sum <- sum + inv1.list[[i]] %*% cov12.list[[i]] %*% inv2.list[[i]]
}
cov.EB <- sumInv1.inv %*% sum %*% sumInv2.inv
list(beta=beta.EB, cov=cov.EB)
} # END: meta.EB
# Function to get the parms/cov for a study given a named vec
getBetaObject <- function(vec, vars) {
# vec Named vector
# vars Names to extract (must be ordered G, E, GE)
ret <- vec[vars]
ret
} # END: getBetaObject
# Function to get the cov for a study given a named vec
getCovObject <- function(vec, vars) {
# vec Named vector
# vars Names to extract (must be ordered, upper triangle of matrix)
len <- length(vars)
n <- (-1 + sqrt(8*len + 1))/2
ret <- matrix(data=NA, nrow=n, ncol=n)
a <- 1
for (i in 1:n) {
b <- a + n - i
ret[i, i:n] <- vec[vars[a:b]]
ret[i:n, i] <- ret[i, i:n]
a <- b + 1
}
ret
} # END: getCovObject
# Function to get the UML-CML cov for a study given a named vec
getCov12Object <- function(vec, vars) {
# vec Named vector
# vars Names to extract (must be ordered, rows of matrix)
n <- sqrt(length(vars))
ret <- matrix(data=vec[vars], byrow=TRUE, nrow=n, ncol=n)
ret
} # END: getCov12Object
# Function to perfrm meta-analysis for each study
meta.GxE <- function(study.list, op=NULL) {
# study.list List of sublists containing name of file, study name, etc
# op List with names:
# out.file Output file
# snps List of snps to include
# print 0 or 1
##############################################################
# Function to get the number of E parms
##############################################################
get.n.parms <- function(vec, which) {
# vec character vector of column names
# which "E" (for main effects) or "G" (for interactions)
# Remove UML.G.BETA and CML.G.BETA
temp <- !(vec %in% c("UML.G.BETA", "CML.G.BETA"))
vec <- vec[temp]
ids <- grep("BETA", vec, fixed=TRUE)
vec <- vec[ids]
vec <- gsub("UML.", "", vec, fixed=TRUE)
vec <- gsub("CML.", "", vec, fixed=TRUE)
vec <- gsub(".BETA", "", vec, fixed=TRUE)
temp <- substr(vec, 1, 1) == which
vec <- vec[temp]
vec <- unique(vec)
n <- length(vec)
n
} # END: get.n.parms
###############################################################
# Function to get the parm names that would need to be flipped from
# a different reference cat
get.flip.E <- function(vec) {
# Remove UML.G.BETA and CML.G.BETA
temp <- !(vec %in% c("UML.G.BETA", "CML.G.BETA"))
vec <- vec[temp]
ids <- grep("BETA", vec, fixed=TRUE)
vec <- vec[ids]
vec2 <- vec
vec2 <- gsub("UML.", "", vec2, fixed=TRUE)
vec2 <- gsub("CML.", "", vec2, fixed=TRUE)
temp <- substr(vec2, 1, 1) %in% c("E","G")
ret <- vec[temp]
ret
} # END: get.flip.E
###############################################################
# Function to determine which column names would be flipped from
# a different risk allele
get.flip.cnames <- function(vec) {
ids <- grep("BETA", vec, fixed=TRUE)
vec <- vec[ids]
vec2 <- vec
vec2 <- gsub("UML.", "", vec2, fixed=TRUE)
vec2 <- gsub("CML.", "", vec2, fixed=TRUE)
temp <- substr(vec2, 1, 1) == "G"
ret <- vec[temp]
ret
} # END: get.flip.cnames
#################################################################
# Function to flip an allele
flip.allele <- function(a) {
if (a == "A") {
return("T")
} else if (a == "T") {
return("A")
} else if (a == "C") {
return("G")
} else if (a == "G") {
return("C")
} else {
return(NULL)
}
} # END: flip.allele
#################################################################
# Function to determine if G parms need to be flipped
get.flip.status <- function(risk.base, other.base, MAF.base,
risk, other, MAF) {
message <- ""
messageFlag <- 0
flip <- 0
error <- 0
# See if the set of alleles match
set.base <- c(risk.base, other.base)
set <- c(risk, other)
flag <- all(set %in% set.base)
if (!flag) {
# Try flipping set of alleles
set[1] <- flip.allele(set[1])
set[2] <- flip.allele(set[2])
flag <- all(set %in% set.base)
}
if (!flag) {
message <- "Alleles do not match"
messageFlag <- 1
return(flip=0, error=1, message=message)
}
# See if risk allele matches
if (risk != set[1]) flip <- 1
list(flip=flip, error=error, message=message)
} # END: get.flip.status
###################################################################
op <- default.list(op, c("out.file", "print"), list("ERROR", 0),
error=c(1, 0))
nstudy <- length(study.list)
snps <- op[["snps", exact=TRUE]]
snpFlag <- !is.null(snps)
snps.all <- NULL
cnames <- NULL
outfile <- op[["out.file", exact=TRUE]]
outFlag <- !is.null(outfile)
if (snpFlag) snps <- removeWhiteSpace(snps)
# Check each study list
vec <- scan(study.list[[1]]$file, what="character", nlines=1, quiet=TRUE)
for (i in 1:nstudy) {
study.list[[i]] <- check.file.list(study.list[[i]], op=list(vars=vec))
temp <- study.list[[i]][["flip.betas", exact=TRUE]]
if (is.null(temp)) study.list[[i]]$flip.betas <- 0
}
# Read in the data
for (i in 1:nstudy) {
tlist <- study.list[[i]]
x <- loadData.table(tlist)
x[, "SNP"] <- removeWhiteSpace(x[, "SNP"])
if (snpFlag) {
temp <- x[, "SNP"] %in% snps
x <- removeOrKeepRows(x, temp)
}
nx <- nrow(x)
str <- paste("For study ", i, ", ", nx, " SNPs will be included.", sep="")
print(str)
if (!nx) {
str <- paste("Removing study ", i, sep="")
print(str)
study.list[[i]] <- NULL
next
}
rownames(x) <- x[, "SNP"]
tlist$data <- x
study.list[[i]] <- tlist
snps.all <- unique(c(snps.all, x[, "SNP"]))
# Check the names of the columns
if (i == 1) {
cnames <- colnames(x)
} else {
if (!all(cnames == colnames(x))) {
str <- paste("ERROR with column names for study ", i, sep="")
stop(str)
}
}
}
# Get the number of main effect parms and interaction parms
N.E.parms <- get.n.parms(cnames, "E")
N.GE.parms <- get.n.parms(cnames, "G")
# Get the column names for the vector of main effect estimates
UML.beta.cnames <- outNames.me("UML", N.E.parms, N.GE.parms)
CML.beta.cnames <- outNames.me("CML", N.E.parms, N.GE.parms)
UML.cov.cnames <- outNames.cov("UML", N.E.parms, N.GE.parms)
CML.cov.cnames <- outNames.cov("CML", N.E.parms, N.GE.parms)
UML.CML.cov.cnames <- outNames.cov2(N.E.parms, N.GE.parms)
# Get the betas to flip from different reference cat
cat.flip.cnames <- get.flip.E(cnames)
# Get the beta column names to flip for different risk allele
allele.flip.cnames <- get.flip.cnames(cnames)
outnames <- c("SNP", "Risk.allele", "Other.allele", "P.omnibus",
"N.study", "Studies", "Message")
outvec <- character(length(outnames))
names(outvec) <- outnames
# Open the output file
if (outFlag) {
fid <- file(outfile, "w")
writeOut(fid, outnames)
}
# Loop over each snp
nsnps <- length(snps.all)
for (i in 1:nsnps) {
outvec[] <- "NA"
UML.list <- list()
CML.list <- list()
index <- 0
snp <- snps.all[i]
message <- ""
outvec[1:3] <- c(snp, ".", ".")
# For UML-CML cov matrix
UML.CML.list <- list()
# Keep track of alleles and MAF
risk.vec <- NULL
other.vec <- NULL
MAF.vec <- NULL
study.vec <- NULL
# Get the vector of estimates for each study
for (j in 1:nstudy) {
x <- study.list[[j]]$data
if (!(snp %in% rownames(x))) next
vec <- x[snp, ]
names(vec) <- cnames
alleles <- vec["Alleles"]
other <- substr(alleles, 1, 1)
risk <- substr(alleles, 2, 2)
MAF <- as.numeric(vec["MAF"])
vec <- as.numeric(vec[-(1:3)])
if (any(!is.finite(vec))) next
names(vec) <- cnames[-(1:3)]
# Flip betas if needed
if (study.list[[j]]$flip.betas) vec[cat.flip.cnames] <- -vec[cat.flip.cnames]
# Get the UML and CML estimates
CML.beta <- getBetaObject(vec, CML.beta.cnames)
CML.cov <- getCovObject(vec, CML.cov.cnames)
UML.beta <- getBetaObject(vec, UML.beta.cnames)
UML.cov <- getCovObject(vec, UML.cov.cnames)
# For UML-CML cov matrix
UML.CML.cov <- getCov12Object(vec, UML.CML.cov.cnames)
if (any(!is.finite(UML.CML.cov))) next
# Check the alleles
if (index == 0) {
risk.base <- risk
other.base <- other
outvec[2:3] <- c(risk, other)
MAF.base <- MAF
} else {
# Compare to base values
ret <- get.flip.status(risk.base, other.base, MAF.base, risk, other, MAF)
if (ret$error) next
if (ret$flip) {
# Flip betas
UML.beta[allele.flip.cnames] <- -UML.beta[allele.flip.cnames]
CML.beta[allele.flip.cnames] <- -CML.beta[allele.flip.cnames]
}
}
index <- index + 1
#alleles.vec <- c(alleles.vec, alleles)
#MAF.vec <- c(MAF.vec, MAF)
study.vec <- c(study.vec, j)
# Add parms to the lists
tlist <- list(beta=UML.beta, cov=UML.cov)
UML.list[[index]] <- tlist
tlist <- list(beta=CML.beta, cov=CML.cov)
CML.list[[index]] <- tlist
tlist <- list(beta=CML.beta, cov=CML.cov)
UML.CML.list[[index]] <- UML.CML.cov
}
outvec["N.study"] <- index
if (!index) {
outvec["Message"] <- "No studies"
if (outFlag) writeOut(fid, outvec)
next
}
outvec["Studies"] <- paste(study.vec, collapse=",", sep="")
# Perform the fixed effects meta-anlaysis for UML and CML estimates
# Return list: list(meta.beta=meta.beta, meta.cov=meta.cov, meta.cov.inv=sum.sigma.inv,
# inv.list=sigma.inv)
UML.meta <- try(meta.fixed(UML.list), silent=TRUE)
if (checkTryError(UML.meta, conv=0)) {
outvec["Message"] <- "Error in UML meta-analysis"
if (outFlag) writeOut(fid, outvec)
next
}
CML.meta <- try(meta.fixed(CML.list), silent=TRUE)
if (checkTryError(CML.meta, conv=0)) {
outvec["Message"] <- "Error in CML meta-analysis"
if (outFlag) writeOut(fid, outvec)
next
}
# Now perform EB
# meta.EB <- function(beta1, beta2, sumInv1.inv, sumInv2.inv, inv1.list, inv2.list,
# cov12.list)
# Return list: list(beta=beta.EB, cov=cov.EB)
ret <- meta.EB(UML.meta$meta.beta, CML.meta$meta.beta, UML.meta$meta.cov,
CML.meta$meta.cov, UML.meta$inv.list,
CML.meta$inv.list, UML.CML.list)
temp <- waldTest.main(ret$beta, ret$cov, 1:length(UML.beta))
outvec["P.omnibus"] <- temp$pvalue
outvec["Message"] <- "OK"
# Write output
if (outFlag) writeOut(fid, outvec)
}
if (outFlag) close(fid)
outvec
} # END: meta.GxE
Try the CGEN package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.