@@ -1,6 +1,6 @@

 Package: gcatest

 Title: Genotype Conditional Association TEST

-Version: 2.0.2.9000

+Version: 2.0.3.9000

 Author: Wei Hao, Minsun Song, Alejandro Ochoa, John D. Storey

 Maintainer: Alejandro Ochoa <[email protected]>, Wei Hao <[email protected]>, John D. Storey <[email protected]>

 Encoding: UTF-8

@@ -1,5 +1,6 @@

 # Generated by roxygen2: do not edit by hand

+export(delta_deviance_lf)

 export(gcat)

 export(gcat.stat)

 export(gcatest)

@@ -33,4 +33,11 @@ Internally there was major code restructuring, and added unit tests for all func

 - Added internal tests for deviance calculations against `stats::glm`.

 - Deviance code (internal `delta_deviance_snp`) now returns `NA` instead of stopping when an "impossible" case is encountered (when the genotype `x` is non-zero but the fitted probabilities under either null or alternative model are zero, or the alternative allele dosage (`x-2`) has the same problem).

-  These cases are clearly model fitting failures, and can arise for common ill-defined problems, particularly under binary `adjustment` variables passed to `gcat` together with rare variants; these individual cases are not handled any better by `stats:glm`, so it seemed most sensible to return `NA` at such loci and not stop.

+  These cases are clearly model fitting failures, and can arise for common ill-defined problems, particularly under binary `adjustment` variables passed to `gcat` together with rare variants; these individual cases are not handled any better by `stats::glm`, so it seemed most sensible to return `NA` at such loci and not stop.

+

+# 2021-05-11 - gcatest 2.0.3.9000

+

+- Added function `delta_deviance_lf`, which calculates the delta deviance from two logistic models and the genotype matrix data.

+  This function is a more general version of `gcat.stat` (which uses the new function internally), to essentially consider models that differ by more than one degree of freedom.

+  It was written in particular for an external application in mind, namely the `jackstraw` package.

+- Internal function `assoc_snp` was renamed to `delta_deviance_snp_lf` and its last argument changed to match that of `delta_deviance_lf` (alternative logistic factors instead of trait).

new file mode 100644

@@ -0,0 +1,68 @@

+#' Calculate delta deviance of logistic fits of a null and alternative model

+#'

+#' This function fits, at each locus of a given genotype matrix, two logistic models, and under the assumption that the models are nested, calculates the delta deviance between the two.

+#' This general function is intended for testing models in a broad setting; for the specific problem of genetic association, the interface in [gcat()] and [gcat.stat()] are more user-friendly.

+#'

+#' @inheritParams lfa::lfa

+#' @param LF0 Logistic factors for null model.

+#' @param LF1 Logistic factors for alternative model.

+#'

+#' @return The vector of delta deviance values, one per locus of `X`.

+#'

+#' @examples

+#' \dontrun{

+#' devdiff <- delta_deviance_lf( X, LF0, LF1 )

+#' }

+#'

+#' @export

+delta_deviance_lf <- function( X, LF0, LF1 ){

+    if ( missing( X ) )

+        stop( 'Genotype matrix `X` is required!' )

+    if ( missing( LF0 ) )

+        stop( "`LF0` matrix is required!" )

+    if ( missing( LF1 ) )

+        stop( "`LF1` matrix is required!" )

+    

+    # check class

+    is_BEDMatrix <- FALSE

+    if ( "BEDMatrix" %in% class(X) ) {

+        is_BEDMatrix <- TRUE

+    } else if ( !is.matrix( X ) )

+        stop( '`X` must be a matrix!' )

+

+    # get dimensions

+    if ( is_BEDMatrix ) {

+        m <- ncol(X)

+        n <- nrow(X)

+    } else {

+        n <- ncol(X)

+        # m not used in this case

+    }

+

+    # check dimensions

+    if ( nrow(LF0) != n )

+        stop( 'Number of individuals in `X` must equal number of individuals (rows) in `LF0`' )

+    if ( nrow(LF1) != n )

+        stop( 'Number of individuals in `X` must equal number of individuals (rows) in `LF1`' )

+    

+    # check LFs for missing values

+    if ( anyNA( LF0 ) )

+        stop( '`LF0` must not have missing values!' )

+    if ( anyNA( LF1 ) )

+        stop( '`LF1` must not have missing values!' )

+

+    # start processing

+    if ( is_BEDMatrix ) {

+        # explicit loop for BEDMatrix

+        devdiff <- vector('numeric', m)

+        for ( i in 1 : m ) {

+            # get locus i genotype vector

+            xi <- X[ , i ]

+            # calculate and store result

+            devdiff[ i ] <- delta_deviance_snp_lf( xi, LF0, LF1 )

+        }

+    } else

+        devdiff <- apply( X, 1, delta_deviance_snp_lf, LF0, LF1 )

+    return( devdiff )

+}

+

similarity index 55%

rename from R/assoc_snp.R

rename to R/delta_deviance_snp_lf.R

@@ -1,35 +1,35 @@

-assoc_snp <- function(xi, LF, trait){

+delta_deviance_snp_lf <- function( xi, LF0, LF1 ){

     if ( missing( xi ) )

         stop( 'Genotype vector `xi` is required!' )

-    if ( missing( LF ) )

-        stop( "`LF` matrix is required!" )

-    if ( missing( trait ) )

-        stop( "`trait` is required!" )

+    if ( missing( LF0 ) )

+        stop( "`LF0` matrix is required!" )

+    if ( missing( LF1 ) )

+        stop( "`LF1` matrix is required!" )

     # check dimensions

     n <- length( xi )

-    if ( nrow(LF) != n )

-        stop( 'Number of individuals in `xi` must equal number of individuals (rows) in `LF`' )

-    if ( length(trait) != n )

-        stop( 'Number of individuals in `xi` must equal number of individuals in `trait`' )

+    if ( nrow(LF0) != n )

+        stop( 'Number of individuals in `xi` must equal number of individuals (rows) in `LF0`' )

+    if ( nrow(LF1) != n )

+        stop( 'Number of individuals in `xi` must equal number of individuals (rows) in `LF1`' )

     # check trait and LF for missing values

-    if ( anyNA( trait ) )

-        stop( '`trait` must not have missing values!' )

-    if ( anyNA( LF ) )

-        stop( '`LF` must not have missing values!' )

+    if ( anyNA( LF0 ) )

+        stop( '`LF0` must not have missing values!' )

+    if ( anyNA( LF1 ) )

+        stop( '`LF1` must not have missing values!' )

     # remove individuals whose genotypes are NA

     # though `lfa::af_snp` allows NAs, the resulting imputed allele frequencies do not contribute to the deviance (because the xi's are factors of the sum)!  So best to exclude them from the beginning.

     indexes_keep <- !is.na(xi)

     xi <- xi[indexes_keep]

-    LF <- LF[indexes_keep,]

-    trait <- trait[indexes_keep]

-

-    # perform two logistic regressions, under the null and alternative (add trait to LF/covariates), resulting in estimated allele frequencies under each model

-    p0 <- lfa::af_snp( xi, LF )

-    p1 <- lfa::af_snp( xi, cbind(LF, trait) )

-

+    LF0 <- LF0[indexes_keep,]

+    LF1 <- LF1[indexes_keep,]

+    

+    # perform two logistic regressions, under the null and alternative, resulting in estimated allele frequencies under each model

+    p0 <- lfa::af_snp( xi, LF0 )

+    p1 <- lfa::af_snp( xi, LF1 )

+    

     # now compute delta deviance

     # uses a special numerically stable algorithm

     devdiff <- delta_deviance_snp( xi, p0, p1 )

@@ -11,21 +11,13 @@ gcat.stat <- function( X, LF, trait, adjustment = NULL ){

     if ( missing( trait ) )

         stop( "`trait` is required!" )

-    # check class

-    is_BEDMatrix <- FALSE

+    # check class, get dimensions

     if ( "BEDMatrix" %in% class(X) ) {

-        is_BEDMatrix <- TRUE

-    } else if ( !is.matrix( X ) )

-        stop( '`X` must be a matrix!' )

-

-    # get dimensions

-    if ( is_BEDMatrix ) {

-        m <- ncol(X)

         n <- nrow(X)

-    } else {

+    } else if ( is.matrix( X ) ) {

         n <- ncol(X)

-        # m not used in this case

-    }

+    } else 

+        stop( '`X` must be a matrix!' )

     # check dimensions

     if ( nrow(LF) != n )

@@ -48,18 +40,14 @@ gcat.stat <- function( X, LF, trait, adjustment = NULL ){

         # if all good, add adjustments to LFs

         LF <- cbind( LF, adjustment )

     }

+

+    # bind matrices

+    LF1 <- cbind(LF, trait)

+

+    # this performs the calculations in a broader setting

+    # (shared with `jackstraw` package)

+    devdiff <- delta_deviance_lf( X, LF, LF1 )

-    if ( is_BEDMatrix ) {

-        # explicit loop for BEDMatrix

-        devdiff <- vector('numeric', m)

-        for ( i in 1 : m ) {

-            # get locus i genotype vector

-            xi <- X[ , i ]

-            # calculate and store result

-            devdiff[ i ] <- assoc_snp( xi, LF, trait )

-        }

-    } else

-        devdiff <- apply( X, 1, assoc_snp, LF, trait )

     return( devdiff )

 }

new file mode 100644

@@ -0,0 +1,31 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/delta_deviance_lf.R

+\name{delta_deviance_lf}

+\alias{delta_deviance_lf}

+\title{Calculate delta deviance of logistic fits of a null and alternative model}

+\usage{

+delta_deviance_lf(X, LF0, LF1)

+}

+\arguments{

+\item{X}{A matrix of SNP genotypes, i.e. an integer matrix of 0's,

+1's, 2's and \code{NA}s.

+BEDMatrix is supported.

+Sparse matrices of class Matrix are not supported (yet).}

+

+\item{LF0}{Logistic factors for null model.}

+

+\item{LF1}{Logistic factors for alternative model.}

+}

+\value{

+The vector of delta deviance values, one per locus of \code{X}.

+}

+\description{

+This function fits, at each locus of a given genotype matrix, two logistic models, and under the assumption that the models are nested, calculates the delta deviance between the two.

+This general function is intended for testing models in a broad setting; for the specific problem of genetic association, the interface in \code{\link[=gcat]{gcat()}} and \code{\link[=gcat.stat]{gcat.stat()}} are more user-friendly.

+}

+\examples{

+\dontrun{

+devdiff <- delta_deviance_lf( X, LF0, LF1 )

+}

+

+}

@@ -169,11 +169,12 @@ test_that( "delta_deviance_snp works", {

     expect_true( is.na( delta_deviance_snp( xi, pi, pi_bad2 ) ) )

 })

-test_that("assoc_snp works", {

+test_that("delta_deviance_snp_lf works", {

     # test a single SNP

     i <- 1

+    LF1 <- cbind(LFs, trait)

     expect_silent(

-        devdiff <- assoc_snp( X[ i, ] , LFs, trait )

+        devdiff <- delta_deviance_snp_lf( X[ i, ] , LFs, LF1 )

     )

     expect_equal( length( devdiff ), 1 )

     expect_true( is.numeric( devdiff ) )

@@ -188,6 +189,28 @@ test_that("assoc_snp works", {

     expect_equal( devdiff, devdiff_glm )

 })

+test_that("delta_deviance_lf works", {

+    LF1 <- cbind(LFs, trait)

+    expect_silent(

+        devdiff <- delta_deviance_lf( X , LFs, LF1 )

+    )

+    expect_equal( length( devdiff ), m_loci )

+    expect_true( is.numeric( devdiff ) )

+    # delta deviances can be NA if LFA/glm.fit fail to converge

+    expect_true( !any( is.na( devdiff ) ) )

+    # theoretically this is true, but in practice it depends on LFA fitting these models well, so testing this is not appropriate here (fails sometimes)

+    ## expect_true( all( devdiff >= 0 ) )

+    

+    # loop through `glm` version

+    devdiff_glm <- vector( 'numeric', m_loci )

+    for ( i in 1 : m_loci ) {

+        suppressWarnings(

+            devdiff_glm[i] <- delta_deviance_snp_glm( X[ i, ], LFs, trait )

+        )

+    }

+    expect_equal( devdiff, devdiff_glm )

+})

+

 test_that("gcat.stat works", {

     expect_silent(

         devdiff <- gcat.stat( X, LFs, trait )

@@ -253,6 +276,15 @@ if (

     # load as a BEDMatrix object

     X_BEDMatrix <- suppressMessages(suppressWarnings( BEDMatrix( file_bed ) ))

+    test_that("delta_deviance_lf works with BEDMatrix", {

+        LF1 <- cbind( LFs, trait )

+        devdiff_basic <- delta_deviance_lf( X, LFs, LF1 )

+        expect_silent(

+            devdiff_BM <- delta_deviance_lf( X_BEDMatrix, LFs, LF1 )

+        )

+        expect_equal( devdiff_basic, devdiff_BM )

+    })

+

     test_that("gcat.stat works with BEDMatrix", {

         devdiff_basic <- gcat.stat( X, LFs, trait )

         expect_silent(