[Gtdb-commits] r52 - in pkg/gt.db: R man

[noreply at r-forge.r-project.org]

Author: dahinds
Date: 2010-05-24 22:22:08 +0200 (Mon, 24 May 2010)
New Revision: 52

Modified:
   pkg/gt.db/R/score.R
   pkg/gt.db/man/score.glm.Rd
Log:
- Added choice of scoring functions (LR, Wald, Rao) to score.glm



Modified: pkg/gt.db/R/score.R
===================================================================
--- pkg/gt.db/R/score.R	2010-05-20 20:50:07 UTC (rev 51)
+++ pkg/gt.db/R/score.R	2010-05-24 20:22:08 UTC (rev 52)
@@ -84,7 +84,7 @@
     )), fit='jt')
 }
 
-score.chisq.2x2 <- function(formula, data, ploidy)
+score.chisq.2x2 <- function(formula, data, ploidy=NA)
 {
     pt <- data[,.lhs.formula(formula)]
     gt <- factor(data$genotype,levels=0:2)
@@ -227,14 +227,23 @@
 #-----------------------------------------------------------------------
 
 score.glm <-
-function(formula, data, ploidy, drop='genotype',
+function(formula, data, ploidy, drop='genotype', test=c('LR','Wald','Rao'),
          mode=c('additive','recessive','dominant','general'))
 {
     data$genotype <- .recode.gt(data$genotype, match.arg(mode))
     m <- glm(formula, binomial, data)
-    nf <- .null.model(formula, term=drop)
-    pvalue <- anova(update(m,nf), m, test='Chisq')[2,5]
     x <- coef(summary(m))
+    test <- match.arg(test)
+    if (test == 'LR') {
+        nf <- .null.model(formula, term=drop)
+        pvalue <- anova(update(m,nf), m, test='Chisq')[2,5]
+    } else if (test == 'Rao') {
+        require(statmod)
+        nf <- .null.model(formula, term=drop)
+        pvalue <- 2*pnorm(-abs(glm.scoretest(update(m,nf), data[,drop])))
+    } else {
+        pvalue <- x[drop,4]
+    }
     x <- if (drop %in% row.names(x)) unname(x[drop,]) else c(NA,NA)
     keep.attr(data.frame(
         pvalue, effect=x[1], stderr=x[2]

Modified: pkg/gt.db/man/score.glm.Rd
===================================================================
--- pkg/gt.db/man/score.glm.Rd	2010-05-20 20:50:07 UTC (rev 51)
+++ pkg/gt.db/man/score.glm.Rd	2010-05-24 20:22:08 UTC (rev 52)
@@ -24,6 +24,7 @@
 }
 \usage{
 score.glm(formula, data, ploidy, drop='genotype',
+          test=c('LR','Wald','Rao'),
           mode=c('additive','recessive','dominant','general'))
 }
 \arguments{
@@ -31,13 +32,18 @@
   \item{data}{a data frame containing the variables in the model.}
   \item{ploidy}{see \code{\link{fetch.gt.data}}.  Ignored.}
   \item{drop}{the model term to test for association.}
+  \item{test}{the scoring test.  See details.}
   \item{mode}{genetic mode of action to be tested.}
 }
 \details{
   The model formula may include covariates and should contain a single
-  genotype term without interactions.  Significance is assessed by
-  ANOVA comparing the full model with a null model constructed by
-  removing this term. 
+  genotype term without interactions.  For \code{test='LR'} (the default),
+  significance is assessed by ANOVA comparing the full model with a null
+  model constructed by removing this term.  Other options are
+  \code{test='Wald'} (a Wald test on the coefficient of the genotype
+  term in the regression model), and \code{test='Rao'} (the Rao score
+  test, which is equivalent to the Cochran-Armitage trend test if there
+  are no covariates).
 
   At sex linked loci, haploid males are treated the same as the
   corresponding diploid female homozygotes.