[GenABEL-dev] Patch bug [#1322] polygenic would not run with a covariate which is non-informative

[pirastu at burlo.trieste.it]

Dear all,

attached you can find a patch to fix this bug:
[#1322] polygenic would not run with a covariate which is non-informative

Basically if any of the covariates is non informative it droppes it and
gives a warning of which variables it has dropped. If non covariates are
left it will stop and warn you about that too. The output gives also
another warning but I think this is dependant of the new version of the
function and not by my changes.
Here is the script I've tested with the output:

library(GenABEL)
data(ge03d2.clean)
df <- ge03d2.clean[1:200,autosomal(ge03d2.clean)]
gkin <- ibs(df,w="freq")
# get a SNP which is informative + no missings
s <- summary(gtdata(df))
maf <- pmin(s[,"Q.2"],1.-s[,"Q.2"])
nomisSnp <- (rownames(s)[which(maf>0.2 & s[,"CallRate"]==1)])[1]
snp <- as.numeric(gtdata(df[, nomisSnp]))
table(snp,exclude=NULL)
pol1 <- polygenic(height~sex+age+snp,df,kin=gkin,quiet=TRUE)
pol1$esth2
# introduce missing indicator
snpInd <- rep(1,length(snp))
snpInd[is.na(snp)] <- NA
# this does did not work!

pol0 <- polygenic(height~sex+age+snpInd,df,kin=gkin,quiet=TRUE)

#now it does
Warning in polygenic(height ~ sex + age + snpInd, df, kin = gkin, quiet =
TRUE) :
  some eigenvalues <=0, taking ABS for det; try option llfun='polylik'
Warning message:
In polygenic(height ~ sex + age + snpInd, df, kin = gkin, quiet = TRUE) :
  Variable snpInd has only one value: Dropped

str(pol0)
List of 8
 $ h2an       :List of 5
  ..$ minimum   : num 1012
  ..$ estimate  : num [1:5] 170.691 12.537 -0.172 0.122 58.121
  ..$ gradient  : num [1:5] -2.70e-09 -3.79e-09 -2.61e-07 -4.97e-10 -1.21e-09
  ..$ code      : int 1
  ..$ iterations: int 57
 $ residualY  : Named num [1:200] -10.01 -3.66 5.47 11.46 -4.92 ...
  ..- attr(*, "names")= chr [1:200] "id4" "id10" "id25" "id33" ...
 $ esth2      : num 0.122
 $ InvSigma   : num [1:200, 1:200] 1.75e-02 -2.09e-05 4.90e-05 -1.22e-05
1.16e-04 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:200] "id4" "id10" "id25" "id33" ...
  .. ..$ : chr [1:200] "id4" "id10" "id25" "id33" ...
 $ measuredIDs: logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
 $ pgresidualY: Named num [1:200] -8.89 -3.25 5.35 10.38 -4.11 ...
  ..- attr(*, "names")= chr [1:200] "id4" "id10" "id25" "id33" ...
 $ call       : language polygenic(formula = height ~ sex + age + snpInd,
kinship.matrix = gkin,      data = df, quiet = TRUE)
 $ convFGLS   : logi TRUE
 - attr(*, "class")= chr "polygenic"

I've also attached a first version of the "How to Make a patch"
intructions. I've based this on what I've actually done to make this patch
and if you have any suggestions comments changes ecc. please let me know.

Best

Nicola

-------------- next part --------------
A non-text attachment was scrubbed...
Name: polygenic_patch
Type: application/octet-stream
Size: 1428 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110524/e55cea80/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: How to patch
Type: application/octet-stream
Size: 2081 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110524/e55cea80/attachment-0001.obj>