assocCoxPH: Cox proportional hazards

Description

Fits Cox proportional hazards model

Usage

assocCoxPH(genoData, event, time.to.event, gene.action = c("additive", "dominant", "recessive"), covar = NULL, ivar = NULL, strata = NULL, scan.exclude = NULL, effectAllele = c("minor", "alleleA"), snpStart = NULL, snpEnd = NULL, block.size = 5000, verbose = TRUE)

Arguments

genoData

a GenotypeData object

event

name of scan annotation variable in genoData for event to analyze (chould be coded 0/1 or FALSE/TRUE)

time.to.event

name of scan annotation variable in genoData for time to event

gene.action

"additive" coding sets the marker variable for homozygous minor allele samples = 2, heterozygous samples = 1, and homozygous major allele samples = 0. "dominant" coding sets the marker variable for homozygous minor allele samples = 2, heterozygous samples = 2, and homozygous major allele samples = 0. "recessive" coding sets the marker variable for homozygous minor allele samples = 2, heterozygous samples = 0, and homozygous major allele samples = 0. (If effectAllele="alleleA", the coding reflects alleleA instead of the minor allele.)

covar

a vector of the names of the covariates to adjust for (columns in the scan annotation of genoData)

ivar

the name of the variable in covar to include as an interaction with genotype

strata

a vector of names of variables to stratify on for a stratified analysis

scan.exclude

a vector of scanIDs for scans to exclude

effectAllele

whether the effects should be returned in terms of the minor allele for the tested sample (effectAllele="minor") or the allele returned by getAlleleA(genoData) (effectAllele="alleleA"). If the minor allele is alleleB for a given SNP, the difference between these two options will be a sign change for the beta estimate.

snpStart

index of the first SNP to analyze, defaults to first SNP

snpEnd

index of the last SNP to analyze, defaults to last SNP

block.size

number of SNPs to read in at once

verbose

logical for whether to print status updates

Value

snpID: the snpIDs
chr: chromosome SNPs are on
n.events: number of events in complete cases for each SNP
effect.allele: which allele ("A" or "B") is the effect allele
EAF: effect allele frequency
MAF: minor allele frequency
filter: TRUE if SNP passes the MAF filter (2*MAF*(1-MAF)*n.events > 75)
Est: beta estimate for genotype
SE: standard error of beta estimate for the genotype
z.Stat: z statistic for association
z.pval: p-value for association
GxE.Stat: Likelihood ratio test statistic for the genotype*ivar interaction parameter
GxE.pval: p-value for the likelihood ratio test statistic

Details

This function performs Cox proportional hazards regression of a survival object (using the Surv function) on SNP genotype and other covariates. It uses the coxph function from the R survival library. It is recommended to filter results returned using 2*MAF*(1-MAF)*n.events > 75 where MAF = minor allele frequency and n.events = number of events. This filter was suggested by Ken Rice and Thomas Lumley, who found that without this requirement, at threshold levels of significance for genome-wide studies, Cox regression p-values based on standard asymptotic approximations can be notably anti-conservative.

Note: Y chromosome SNPs must be analyzed separately because they only use males.

Examples

Run this code

library(GWASdata)
data(illuminaScanADF)
scanAnnot <- illuminaScanADF

# exclude duplicated subjects
scan.exclude <- scanAnnot$scanID[scanAnnot$duplicated]

# create some variables for the scans
scanAnnot$sex <- as.factor(scanAnnot$sex)
scanAnnot$age <- rnorm(nrow(scanAnnot), mean=40, sd=10)
scanAnnot$event <- rbinom(nrow(scanAnnot), 1, 0.4)
scanAnnot$ttoe <- rnorm(nrow(scanAnnot), mean=100, sd=10)

# create data object
gdsfile <- system.file("extdata", "illumina_geno.gds", package="GWASdata")
gds <- GdsGenotypeReader(gdsfile)
genoData <-  GenotypeData(gds, scanAnnot=scanAnnot)

res <- assocCoxPH(genoData,
                  event="event", time.to.event="ttoe",
  		  covar=c("sex", "age"),
  		  scan.exclude=scan.exclude,
  		  snpStart=1, snpEnd=100)

close(genoData)