From d43b2b9c48c1d31cccd6556c452004fd0ade5a48 Mon Sep 17 00:00:00 2001
From: Karl Broman <kbroman@gmail.com>
Date: Sat, 3 Jul 2021 07:51:08 -0500
Subject: [PATCH] Add addcovar argument to cim()

 - Implements Issue #95
---
 DESCRIPTION |  4 ++--
 NEWS.md     |  8 ++++++++
 R/cim.R     | 22 +++++++++++++++-------
 man/cim.Rd  |  3 ++-
 4 files changed, 27 insertions(+), 10 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 0dfb063ed..37fa2d549 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: qtl
-Version: 1.48-1
-Date: 2021-03-24
+Version: 1.49-1
+Date: 2021-06-30
 Title: Tools for Analyzing QTL Experiments
 Author: Karl W Broman <broman@wisc.edu> and Hao Wu, with
   ideas from Gary Churchill and Saunak Sen and contributions from
diff --git a/NEWS.md b/NEWS.md
index 923750f8f..787ceaa97 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -16,6 +16,14 @@ copyright (c) 2001-2021, Karl W Broman
     at https://www.r-project.org/Licenses/GPL-3
 ----------------------------------------------------------------------
 
+## Verison 1.49, 2021-06-30
+
+## New features
+
+- `cim()` now includes an `addcovar` argument for including additional
+  covariates in the model.
+
+
 ## Version 1.48, 2021-03-24
 
 ### Bug fixes
diff --git a/R/cim.R b/R/cim.R
index 02bfe2211..88d33830c 100644
--- a/R/cim.R
+++ b/R/cim.R
@@ -2,8 +2,8 @@
 #
 # cim.R
 #
-# copyright (c) 2007-2019, Karl W Broman
-# last modified Dec, 2019
+# copyright (c) 2007-2021, Karl W Broman
+# last modified Jun, 2021
 # first written Jan, 2007
 #
 #     This program is free software; you can redistribute it and/or
@@ -34,12 +34,19 @@ cim <-
              method=c("em", "imp", "hk", "ehk"),
              imp.method=c("imp", "argmax"), error.prob=0.0001,
              map.function=c("haldane", "kosambi", "c-v", "morgan"),
-             n.perm)
+             addcovar=NULL, n.perm)
 {
     method <- match.arg(method)
     imp.method <- match.arg(imp.method)
     map.function <- match.arg(map.function)
 
+    if(!is.null(addcovar)) {
+        if(!is.matrix(addcovar)) addcovar <- as.matrix(addcovar)
+        if(nrow(addcovar) != nind(cross)) {
+            stop("nrow(addcovar) != no. individuals in cross")
+        }
+    }
+
     type <- crosstype(cross)
     if(type=="4way")
         stop("cim() has not been implemented for 4-way crosses")
@@ -78,7 +85,8 @@ cim <-
             cross$pheno[,pheno.col] <- y
             temp <- cim(cross, pheno.col=pheno.col, n.marcovar=n.marcovar,
                         window=window, method=method, imp.method=imp.method,
-                        error.prob=error.prob, map.function=map.function)
+                        error.prob=error.prob, map.function=map.function,
+                        addcovar=addcovar)
             results[i,1] <- max(temp[,3], na.rm=TRUE)
         }
         class(results) <- c("scanoneperm", "matrix")
@@ -104,7 +112,7 @@ cim <-
     if(type=="f2") useac <- expandf2covar(ac)
     else useac <- ac
 
-    firstscan <- scanone(cross, pheno.col=pheno.col, addcovar=useac,
+    firstscan <- scanone(cross, pheno.col=pheno.col, addcovar=cbind(useac,addcovar),
                          method=method)
 
     # scan again, dropping one marker covariate at a time
@@ -112,7 +120,7 @@ cim <-
         if(type=="f2") useac <- expandf2covar(ac[,-i])
         else useac <- ac[,-i]
 
-        temp <- scanone(cross, pheno.col=pheno.col, addcovar=useac,
+        temp <- scanone(cross, pheno.col=pheno.col, addcovar=cbind(useac,addcovar),
                         method=method, chr=chrpos[i,1])
         wh1 <- (firstscan[,1]==chrpos[i,1] & firstscan[,2] >= chrpos[i,2]-window &
                 firstscan[,2] <= chrpos[i,2]+window)
@@ -149,7 +157,7 @@ cim <-
 
                     if(type=="f2") useac <- expandf2covar(tempac)
                     else useac <- tempac
-                    temp <- scanone(cross, pheno.col=pheno.col, addcovar=useac,
+                    temp <- scanone(cross, pheno.col=pheno.col, addcovar=cbind(useac,addcovar),
                                     method=method, chr=j)
 
                     firstscan[whpos2,3] <- temp[whpos,3]
diff --git a/man/cim.Rd b/man/cim.Rd
index 62e584e17..5a4a1cdaf 100644
--- a/man/cim.Rd
+++ b/man/cim.Rd
@@ -16,7 +16,7 @@ cim(cross, pheno.col=1, n.marcovar=3, window=10,
     method=c("em", "imp", "hk", "ehk"),
     imp.method=c("imp", "argmax"), error.prob=0.0001,
     map.function=c("haldane", "kosambi", "c-v", "morgan"),
-    n.perm)
+    addcovar=NULL, n.perm)
 }
 \arguments{
  \item{cross}{An object of class \code{cross}. See
@@ -39,6 +39,7 @@ cim(cross, pheno.col=1, n.marcovar=3, window=10,
    the missing marker genotype data.}
  \item{map.function}{Map function used when imputing the missing marker
    genotype data.}
+ \item{addcovar}{Optional numeric matrix of additional covariates to include.}
  \item{n.perm}{If specified, a permutation test is performed rather than
    an analysis of the observed data.  This argument defines the number
    of permutation replicates.}