Missing genotype compatible

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -1 +1,2 @@
 inst/doc
+.Rproj.user

DESCRIPTION

@@ -8,5 +8,9 @@ Maintainer: Diptavo Dutta <diptavo@umich.edu>
 Description: Kernel Regression based association tests for Multiple phenotypes. The functions aggregate variant-phenotype score statistic in a particular region and computes corresponding p-values efficiently.
 Depends: SKAT, nlme, copula
 License: GPL (>=2)
-Suggests: knitr,rmarkdown,R.rsp
+Suggests: 
+    knitr,
+    rmarkdown,
+    R.rsp,
+    testthat
 VignetteBuilder: knitr

MultiSKAT.Rproj

@@ -0,0 +1,17 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

R/Genotype.Kernels.R

@@ -11,9 +11,9 @@ function(Z,obj.res, kernel = "linear.weighted", Is.Common=FALSE, weights.beta=c(
   if(m < 3 & verbose){
 	msg <-sprintf("Rare variant test with < 3 variants is not advisable")
 	warning(msg,call. = FALSE)}
-  n.rare <- sum(MAF(Z) < 0.01)
+  n.rare <- sum(MAF(Z) < 0.01, na.rm = TRUE)
     if(verbose){print(paste("The region has ",n.rare," rare variants",sep = ""))}
-  mc <- sum(MAC(Z))
+  mc <- sum(MAC(Z), na.rm = TRUE)
   if(mc < 5 & verbose){
 	msg <-sprintf("Rare variant test with total MAC < 5 is not advisable")
 	warning(msg,call. = FALSE)}

R/Helper.R

@@ -1,28 +1,28 @@
-mat.sqrt <- function(A)
-{
-  ei<-eigen(A)
-  d<-ei$values
-  d<-(d+abs(d))/2
-  d2<-sqrt(d)
-  ans<-ei$vectors %*% diag(d2) %*% t(ei$vectors)
-  return(ans)
-}
-
-MAF <- function(G){
-  mf <- array()
-  for(i in 1:ncol(G)){
-    mf[i] <- sqrt(sum(G[,i] == 0)/nrow(G));
-    if(mf[i] > 0.5) 
-      mf[i] = 1- mf[i];
-  }
-  return(mf)
-}
-
-MAC <- function(G){
-  mc <- array()
-  for(i in 1:ncol(G)){
-    mc[i] <- sum(G[,i] == 1) + 2*min(sum(G[,i] == 2),sum(G[,i] == 0));
-  }
-  return(mc);
-}
-
+mat.sqrt <- function(A)
+{
+  ei<-eigen(A)
+  d<-ei$values
+  d<-(d+abs(d))/2
+  d2<-sqrt(d)
+  ans<-ei$vectors %*% diag(d2) %*% t(ei$vectors)
+  return(ans)
+}
+
+MAF <- function(G){
+  mf <- array()
+  for(i in 1:ncol(G)){
+    mf[i] <- sqrt(sum(G[,i] == 0, na.rm = TRUE)/nrow(G));
+    if(mf[i] > 0.5) 
+      mf[i] = 1- mf[i];
+  }
+  return(mf)
+}
+
+MAC <- function(G){
+  mc <- array()
+  for(i in 1:ncol(G)){
+    mc[i] <- sum(G[,i] == 1, na.rm = TRUE) + 2*min(sum(G[,i] == 2, na.rm = TRUE), sum(G[,i] == 0, na.rm = TRUE), na.rm = TRUE);
+  }
+  return(mc);
+}
+

tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(MultiSKAT)
+
+test_check("MultiSKAT")

tests/testthat/test-missing-genotype-data.R

@@ -0,0 +1,18 @@
+context("test-missing-genotype-data")
+
+
+
+test_that("multiplication works", {
+
+  data(MultiSKAT.example.data)
+  attach(MultiSKAT.example.data)
+
+  ## Introduce NAs in Genotypes
+  Genotypes[1,1] <- NA
+
+  expect_warning(MultiSKAT(MultiSKAT_NULL(Phenotypes, Cov), Genotypes, Sigma_p = cov(Phenotypes), verbose = FALSE), 
+                 "The missing genotype rate is 0.000004. Imputation is applied.")
+
+  detach(MultiSKAT.example.data)
+
+})