[GenABEL-dev] faster polygenic

[Yurii Aulchenko]

Hi Xia,

just gave your procedure a try -- and time-wise it works _incredibly_ fast
(code attached FYI)! So, this sounds like a good way to go to speed up
'polygenic' -- I thinking of writing a wrapper (say, 'polygenic_hglm') which
will require(hglm), prepare the data, call hglm and the format the output to
make a valid 'polygenic' object.

Few things I could not see answer right away -- can you please help? -- it
will save me quite some time!

1) do I understand correctly that fixed effects estimates are in
hglm(...)$fixef
2) do I understand correctly that 'breeding values' are in hglm(...)$ranef
3) do I understand correctly that environmental residuals are in
hglm(...)$resid
4) how can I get hold on the estimates of polygenic and environmental
(residual) variances, and find out what was heritability estimate
5) is there an easy way to get hold on inverse of variance-covariance
matrix. I mean given parameters estimates, it is no problem to compute that,
but just in case it is already hanging around would be stupid to do so

best, and many thanks for you efforts,
Yurii

[code]
library(GenABEL)
library(hglm)
data(ge03d2.clean)
df <- ge03d2.clean[1:200,autosomal(ge03d2.clean)]
gkin <- ibs(df,w="freq")

formula <- as.formula(height ~ sex + age)
t0 <- system.time(res0 <- polygenic(formula,kin=gkin,df))

grel <- gkin
grel[upper.tri(grel)] <- t(grel)[upper.tri(grel)]
data <- phdata(df)
mf <- model.frame(formula,data,na.action=na.omit,drop.unused.levels=TRUE)
y <- model.response(mf)
desmat <- model.matrix(formula,mf)
phids <- rownames(data)[rownames(data) %in% rownames(mf)]
relmat <- grel[phids,phids]*2.0
s <- svd(grel)
L <- s$u %*% diag(sqrt(s$d))
t1 <- system.time(res1 <- hglm(y = y, X = desmat, Z = L, family =
gaussian()))

[/code]

On Wed, Feb 16, 2011 at 2:11 PM, Xia Shen <xia.shen at lcb.uu.se> wrote:

> Hi Yurii,
>
> Bad news and good news:
>
> [BAD]
> I cannot open the two links that you sent...
>
> [GOOD]
> One of my colleague happens to be playing with your 'polygenic' function,
> and has got stuck by some binary trait because the numerical algorithm
> didn't converge. However, treating the trait as normal works. To fit the
> same mixed model, I tried my 'hglm' package (published on the recent issue
> of The R Journal), it worked pretty well, and fitting binary, poisson,
> gamma, ... traits will all work.
>
> The way to fit hglm using GenABEL results is straightforward:
>
> 1. create the IBS matrix G, has to be symmetric, just need to replicate the
> lower triangle of the IBS from GenABEL and overwrite the upper triangle
> 2. create L that satisfies LL' = G by
>    s <- svd(G)
>    L <- s$u%*%diag(sqrt(s$d))
> 3. call hglm like
>    hglm(y = phdata(...)$..., X = design matrix of fixed effects, Z = L,
> family = gaussian())
>    where NA has to be removed before calling, gaussian() can be changed to
> binomial(link = 'logit'), etc.
>
> The hglm function will return all the parameter estimates (fixed effects,
> variance components) and random effects predicts and residuals, etc. The
> algorithm is based on inter-connected GLMs (h-likelihood theory, the work by
> Lee & Nelder in the last decade), efficient and stable, without establishing
> the likelihood and calling numeric algorithm like FGLS.
>
> I think the residuals from the mixed model is the most interesting stuff
> that will be passed into the next GWAS in GenABEL.
>
> Hopefully this can serve some help :)
>
> Cheers
>
> Xia
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110224/debfbc45/attachment.htm>