##########################
## pbat.m(...) examples ##
##########################
## Note, when you run the example (or anything else) you will generally
## get a warning message that the column headers were guessed.
## This means they were guessed, and while I've tried to catch most
## cases, the warning stands for ones I might have missed.
## These cannot be run verbatim, and are just meant to be examples.
##############################
## Further formula examples ##
##############################
# load in the data
# Here we assume that:
# data.phe contains 'preds1', 'preds2', 'preds3', 'time',
# 'censor', 'phenos1', ... 'phenos4'
# data.ped contains 'snp1', 'snp2', 'snp3',
# 'block1snp1','block1snp2',
# 'block2snp1','block2snp2'
data.phe <- read.phe( "data" )
data.ped <- read.ped( "data" )
# This model does just the affection status (always given as
# AffectionStatus) as the phenotype, no predictor covariates, and all
# the snps for a snps analysis.
# Since affection status is dichotomous, we additionally set
# distribution='categorical'
# offset='none'
# NONE is a special keyword to indicate none, and can be only used in
# this case (note that it is _case_ _sensative_);
# otherwise one specifies values from the phenotype object, after and
# including AffectionStatus.
res <- pbat.m( AffectionStatus ~ NONE, phe, ped, fbat="gee",
distribution='categorical', offset='none', ... )
summary( res )
res # equivalent to print(res)
# basic model with one phenotype, does all snps (if none specified)
pbat.m( phenos1 ~ preds1, phe, ped, fbat="gee" )
# same model, but with more phenotypes; here we test them all at once
pbat.m( phenos1 + phenos2 + phenos3 ~ preds1, phe, ped, fbat="gee" )
# same model as just before, but now supposing that these phenotypes are
# instead from a longitudinal study
pbat.m( phenos1 + phenos2 + phenos3 ~ preds1, phe, ped, fbat="pc" )
# like our second model, but the mi() tells it should be a marker
# interaction
pbat.m( phenos1 ~ mi(preds1), phe, ped, fbat="gee" )
# logrank analysis - fbat need not be set
# uses more than one predictor variable
res <- pbat.m( time & censor ~ preds1 + preds2 + preds3, phe, ped )
plot( res )
# single snp analysis (because each snp is seperated by a vertical bar
# '|'), and stratified by group (presence of censor auto-indicates
# log-rank analysis). Note that the group is at the end of the
# expression, and _must_ be at the end of the expression
res <- pbat.m( time & censor ~ preds1^3 + preds2 | snp1 | snp2 |
snp3 / group, temp )
plot( res )
# haplotype analysis, stratified by group
res <- pbat.m( time & censor ~ preds1^2 + preds2^3 | block1snp1
+ block1snp2 | block2snp1 + block2snp2 / group, temp )
# set any of the various options
res <- pbat.m( phenos ~ preds, phe, ped, fbat="pc",
null="linkage, no association", alpha=0.1 )
## New multimarker test (as described above)
# mmaxphenos and mmaxsnps are set to the minimum if not specified
res <- pbat.m( phenos1 + phenos2 + phenos3 ~ preds | m1 | m2 | m3 | m4,
phe, ped, fbat="pc", mminphenos=2, mminsnps=2 )
## And the top markers by conditional power
top( res )
############################
## pbat.obj(...) examples ##
############################
# These will not function; they only serve as examples.
# ... just indicates there are various options to be put here!
res <- pbat.obj("pedfile", snps=c("snp1,snp2"), preds="pred1", ... )
summary(res)
res
# plot is only available for "logrank"
res <- pbat.obj(..., fbat="logrank")
plot( res )
##############################
## pbat.files(...) examples ##
##############################
# These will not function, but only serve as examples.
# Note in the following example, both "pedfile.ped" and "pedfile.phe"
# must exist. If the names differed, then you must specify the
# option 'phe="phefile.phe"', for example.
res <- pbat.files( "pedfile", phenos=c("phenos1","phenos2"),
screening="conditional power" )
summary(res)
resRun the code above in your browser using DataLab