R: Genome-wide FGLS...

GWFGLS {MixABEL}R Documentation

Genome-wide FGLS...

Description

Genome-wide FGLS

Usage

GWFGLS(formula, data, subset, weights, na.action, contrasts, offset, W,
    inverse=TRUE, na.SNP="impute", mincall=0.95, residuals=FALSE,
    test="wald", model.SNP="additive", genodata, gtcoding="typed",
    verbosity=1, varcov=FALSE, include.means=TRUE, singular="ignore",
    with.lm=FALSE, old=FALSE)

Arguments

formula analysis formula; should contain special 'SNP' term
data phenotypic data frame, or 'gwaa.data-class' object
subset subset of data (individuals)
weights RESERVED FOR FURTURE USE
na.action RESERVED FOR FURTURE USE
contrasts RESERVED FOR FURTURE USE
offset RESERVED FOR FURTURE USE
W for GLS, (inverse of) variance-covariance matrix, as such returned by GenABEL's polygenic(...)\$InvSigma, or NULL for LS
inverse wheather W is already inverted
na.SNP how to deal with missing SNP data; 'impute' – substitute with mean, 'drop' – drop rows with missing genotypes
mincall minimall call rate for a SNP (if less, the SNP is dropped)
residuals use residuals for analysis? (faster, but less precise)
test test to be applied, one of 'wald', 'score' or 'robust'
model.SNP SNP model to apply, one of c("additive","dominantB","recessiveB","overdominant","genotypic")
genodata genotypic data; can be missing if data is of 'gwaa.data-class'. Otherwise can be regular matrix or 'databel' matrix
gtcoding one of c("typed","dose","probability") 'typed' – coded with NA, 0, 1, 2
verbosity what to report; 0 – only test stats, 1 – test stats and estimates concerning 'whichtest' parameters, 2 – test stats and estimates of all parameters
varcov wheather var-cov matrix for estimated parameters is to be reported
include.means weather mean values of predictors are to be reported
singular what to do with linear dependencies in X (now only 'ignore' implemented)
with.lm whether LM should be run along (only test purposes; always falls back to 'old' R-only implementation)
old if TRUE, old R-only code implementation is running (testing purposes, slow)

Details

Genome-wide Feasible GLS

Author(s)

Yurii Aulchenko

Examples

library(MASS)
library(mvtnorm)
library(GenABEL)
data(ge03d2.clean)
NIDS <- 100
df <- ge03d2.clean[1:NIDS,autosomal(ge03d2.clean)]
s <- summary(df@gtdata)
maf <- pmin(s$Q.2,1-s$Q.2)
df <- df[,(maf>0.05)]
gkin <- ibs(df[,autosomal(df)],w="freq")

modelh2 <- 0.8
covars <- c("sex","age")
s2 <- 12^2
betas_cov <- c(170,12,0.01)

X <- as.matrix(cbind(rep(1,NIDS),phdata(ge03d2.clean)[1:NIDS,covars]))
grel <- gkin
grel[upper.tri(grel)] <- t(grel)[upper.tri(grel)]
grel <- 2*grel
Y <- as.vector(rmvnorm(1,mean=(X %*% betas_cov),sigma=s2*(modelh2*grel+diag(dim(grel)[1])*(1-modelh2))))
length(Y)
Y[2] <- NA
df@phdata$Y <- Y

mdl <- Y~age+sex
h2 <- polygenic(mdl,data=df,kin=gkin)
h2$h2an
mdl <- Y~age+sex+SNP

gtOld <- df@gtdata[,1:10]
gtReal <- as.double(gtOld)
gtNew <- as(gtReal,"databel")
aIR <- GWFGLS(mdl,data=phdata(df),genodata = gtReal,verbosity=0) #,model="genotypic")
aIO <- GWFGLS(mdl,data=phdata(df),genodata = gtOld, verbosity=0) #,model="genotypic")
aIN <- GWFGLS(mdl,data=phdata(df),genodata = gtNew, verbosity=0) #,model="genotypic")
aWR <- GWFGLS(mdl,data=phdata(df),W=h2$InvSigma,genodata = gtReal,verbosity=0) #,model="genotypic")
aWO <- GWFGLS(mdl,data=phdata(df),W=h2$InvSigma,genodata = gtOld, verbosity=0) #,model="genotypic")
aWN <- GWFGLS(mdl,data=phdata(df),W=h2$InvSigma,genodata = gtNew, verbosity=0) #,model="genotypic")
aIN
aWN
complex_model <- GWFGLS(Y~age+sex*SNP,data=phdata(df),W=h2$InvSigma,
genodata = gtReal,verbosity=2,varcov=TRUE,include.means=TRUE,model="genotypic")
complex_model

[Package MixABEL version 0.0-9 Index]