tests/testthat/test-scan1coef.R
In rqtl/qtl2: Quantitative Trait Locus Mapping in Experimental Crosses
context("effect scan by scan1coef")
convert_probs2qtl2 <-
function(cross)
{
result <- lapply(cross$geno, function(a) {
pr <- aperm(a$prob, c(1,3,2))
rownames(pr) <- paste(1:nrow(pr))
pr })
class(result) <- c("calc_genoprob", "list")
result
}
# calc estimates via lm(), just one chromosome
eff_via_lm <-
function(probs, pheno, addcovar=NULL, intcovar=NULL, weights=NULL,
se=TRUE)
{
npos <- dim(probs)[2]
nind <- length(pheno)
ngen <- dim(probs)[3]
nadd <- ifelse(is.null(addcovar), 0, ncol(addcovar))
nint <- ifelse(is.null(intcovar), 0, ncol(intcovar))
ncoef <- ngen + nadd + (ngen-1)*nint
result <- matrix(nrow=npos, ncol=ncoef)
SEs <- matrix(nrow=npos, ncol=ncoef)
if(is.null(weights)) wts <- rep(1, nind)
else wts <- weights
if(!is.null(intcovar)) intcovar <- as.matrix(intcovar)
for(i in 1:npos) {
X <- probs[,i,]
if(!is.null(addcovar)) X <- cbind(X, addcovar)
if(!is.null(intcovar)) {
for(j in 1:ncol(intcovar))
X <- cbind(X, probs[,i,-1]*intcovar[,j])
}
lm.out <- lm(pheno ~ -1 + X, weights=weights)
result[i,] <- lm.out$coef
if(se) # need to deal with NAs
SEs[i,!is.na(lm.out$coef)] <- summary(lm.out)$coef[,2]
else SEs <- NULL
}
attr(result, "SE") <- SEs
class(result) <- c("scan1coef", "scan1", "matrix")
result
}
test_that("scan1coef for backcross", {
set.seed(9308594)
library(qtl)
data(hyper)
hyper <- hyper[4,] # chr 4 only
hyper <- calc.genoprob(hyper, step=5)
prob <- hyper$geno[[1]]$prob
phe <- hyper$pheno[,1]
weights <- runif(nind(hyper), 0, 5)
covar <- cbind(sex=sample(0:1, nind(hyper), replace=TRUE))
rownames(prob) <- names(phe) <- names(weights) <- rownames(covar) <- paste(1:nind(hyper))
# probs for qtl2 code
prob2 <- convert_probs2qtl2(hyper)
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
coef <- scan1coef(prob2, phe, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA")))
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
coef <- scan1coef(prob2, phe, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
expect_equivalent(coef, expected)
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(nind(hyper)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
})
test_that("scan1coef for backcross, with contrasts", {
skip_on_cran()
set.seed(9308594)
library(qtl)
data(hyper)
hyper <- hyper[4,] # chr 4 only
hyper <- calc.genoprob(hyper, step=5)
prob <- hyper$geno[[1]]$prob
phe <- hyper$pheno[,1]
weights <- runif(nind(hyper), 0, 5)
covar <- cbind(sex=sample(0:1, nind(hyper), replace=TRUE))
rownames(prob) <- names(phe) <- names(weights) <- rownames(covar) <- paste(1:nind(hyper))
# probs for qtl2 code
prob2 <- convert_probs2qtl2(hyper)
# use contrasts
contrasts <- cbind(mu=c(1,1), a=c(-0.5, 0.5))
# change to contrasts for use with lm()
p <- prob
dim(p) <- c(prod(dim(prob)[1:2]), dim(prob)[3])
p <- p %*% contrasts
dim(p) <- dim(prob)
prob <- p
dimnames(prob) <- list(rownames(prob2), dimnames(prob2[[1]])[[3]], colnames(contrasts))
posnames <- dimnames(prob2[[1]])[[3]]
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
coef <- scan1coef(prob2, phe, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a")))
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex", "a:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex", "a:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(nind(hyper)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
})
test_that("scan1coef for intercross", {
skip_on_cran()
set.seed(9308594)
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
prob2 <- calc_genoprob(iron[,7], map, error_prob=0.002) # just look at chr 7
phe <- iron$pheno[,1] # liver phenotype
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
covar <- cbind(sex=covar)
weights <- runif(n_ind(iron), 0, 5)
names(weights) <- names(phe)
# different organization of probs for qtl2 and lm() code
prob <- aperm(prob2[[1]], c(1,3,2)) # rearrange as expected for lm()
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, zerosum=FALSE)
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
posnames <- dimnames(prob2[[1]])[[3]]
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB")))
coef <- scan1coef(prob2, phe, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(n_ind(iron)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
})
test_that("scan1coef for intercross, with contrasts", {
skip_on_cran()
set.seed(9308594)
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
prob2 <- calc_genoprob(iron[,7], map, error_prob=0.002) # just look at chr 7
phe <- iron$pheno[,1] # liver phenotype
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
covar <- cbind(sex=covar)
weights <- runif(n_ind(iron), 0, 5)
names(weights) <- names(phe)
# different organization of probs for qtl2 and lm() code
prob <- aperm(prob2[[1]], c(1,3,2)) # rearrange as expected for lm()
# use contrasts
contrasts <- cbind(mu=c(1,1,1), a=c(-1, 0, 1), d=c(0, 1,0))
# change to contrasts for use with lm()
dn <- dimnames(prob)
p <- prob
dim(p) <- c(prod(dim(prob)[1:2]), dim(prob)[3])
p <- p %*% contrasts
dim(p) <- dim(prob)
prob <- p
dimnames(prob) <- list(dn[[1]], dn[[2]], colnames(contrasts))
posnames <- dimnames(prob2[[1]])[[3]]
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
coef <- scan1coef(prob2, phe, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d")))
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(n_ind(iron)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames,
c("mu", "a", "d", "sex", "another", "a:sex", "d:sex",
"a:another", "d:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
})
test_that("scan1coef deals with mismatching individuals", {
skip_on_cran()
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
probs <- probs[,"3"]
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
phe <- iron$pheno[,2,drop=FALSE]
ind <- c(1:50, 101:150)
expected <- scan1coef(probs[ind,], phe[ind,,drop=FALSE], addcovar=X[ind])
expect_equal(scan1coef(probs[ind,], phe, addcovar=X), expected)
expect_equal(scan1coef(probs, phe[ind,,drop=FALSE], addcovar=X), expected)
expect_equal(scan1coef(probs, phe, addcovar=X[ind]), expected)
})
rqtl/qtl2 documentation built on March 20, 2024, 6:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
context("effect scan by scan1coef")
convert_probs2qtl2 <-
function(cross)
{
result <- lapply(cross$geno, function(a) {
pr <- aperm(a$prob, c(1,3,2))
rownames(pr) <- paste(1:nrow(pr))
pr })
class(result) <- c("calc_genoprob", "list")
result
}
# calc estimates via lm(), just one chromosome
eff_via_lm <-
function(probs, pheno, addcovar=NULL, intcovar=NULL, weights=NULL,
se=TRUE)
{
npos <- dim(probs)[2]
nind <- length(pheno)
ngen <- dim(probs)[3]
nadd <- ifelse(is.null(addcovar), 0, ncol(addcovar))
nint <- ifelse(is.null(intcovar), 0, ncol(intcovar))
ncoef <- ngen + nadd + (ngen-1)*nint
result <- matrix(nrow=npos, ncol=ncoef)
SEs <- matrix(nrow=npos, ncol=ncoef)
if(is.null(weights)) wts <- rep(1, nind)
else wts <- weights
if(!is.null(intcovar)) intcovar <- as.matrix(intcovar)
for(i in 1:npos) {
X <- probs[,i,]
if(!is.null(addcovar)) X <- cbind(X, addcovar)
if(!is.null(intcovar)) {
for(j in 1:ncol(intcovar))
X <- cbind(X, probs[,i,-1]*intcovar[,j])
}
lm.out <- lm(pheno ~ -1 + X, weights=weights)
result[i,] <- lm.out$coef
if(se) # need to deal with NAs
SEs[i,!is.na(lm.out$coef)] <- summary(lm.out)$coef[,2]
else SEs <- NULL
}
attr(result, "SE") <- SEs
class(result) <- c("scan1coef", "scan1", "matrix")
result
}
test_that("scan1coef for backcross", {
set.seed(9308594)
library(qtl)
data(hyper)
hyper <- hyper[4,] # chr 4 only
hyper <- calc.genoprob(hyper, step=5)
prob <- hyper$geno[[1]]$prob
phe <- hyper$pheno[,1]
weights <- runif(nind(hyper), 0, 5)
covar <- cbind(sex=sample(0:1, nind(hyper), replace=TRUE))
rownames(prob) <- names(phe) <- names(weights) <- rownames(covar) <- paste(1:nind(hyper))
# probs for qtl2 code
prob2 <- convert_probs2qtl2(hyper)
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
coef <- scan1coef(prob2, phe, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA")))
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
coef <- scan1coef(prob2, phe, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
expect_equivalent(coef, expected)
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(colnames(prob), c("BB", "BA", "sex", "BA:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(nind(hyper)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(colnames(prob), c("BB", "BA", "sex", "another", "BA:sex", "BA:another")))
})
test_that("scan1coef for backcross, with contrasts", {
skip_on_cran()
set.seed(9308594)
library(qtl)
data(hyper)
hyper <- hyper[4,] # chr 4 only
hyper <- calc.genoprob(hyper, step=5)
prob <- hyper$geno[[1]]$prob
phe <- hyper$pheno[,1]
weights <- runif(nind(hyper), 0, 5)
covar <- cbind(sex=sample(0:1, nind(hyper), replace=TRUE))
rownames(prob) <- names(phe) <- names(weights) <- rownames(covar) <- paste(1:nind(hyper))
# probs for qtl2 code
prob2 <- convert_probs2qtl2(hyper)
# use contrasts
contrasts <- cbind(mu=c(1,1), a=c(-0.5, 0.5))
# change to contrasts for use with lm()
p <- prob
dim(p) <- c(prod(dim(prob)[1:2]), dim(prob)[3])
p <- p %*% contrasts
dim(p) <- dim(prob)
prob <- p
dimnames(prob) <- list(rownames(prob2), dimnames(prob2[[1]])[[3]], colnames(contrasts))
posnames <- dimnames(prob2[[1]])[[3]]
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
coef <- scan1coef(prob2, phe, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a")))
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex", "a:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "a:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "sex", "a:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(nind(hyper)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("mu", "a", "sex", "another", "a:sex", "a:another")))
})
test_that("scan1coef for intercross", {
skip_on_cran()
set.seed(9308594)
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
prob2 <- calc_genoprob(iron[,7], map, error_prob=0.002) # just look at chr 7
phe <- iron$pheno[,1] # liver phenotype
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
covar <- cbind(sex=covar)
weights <- runif(n_ind(iron), 0, 5)
names(weights) <- names(phe)
# different organization of probs for qtl2 and lm() code
prob <- aperm(prob2[[1]], c(1,3,2)) # rearrange as expected for lm()
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, zerosum=FALSE)
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
posnames <- dimnames(prob2[[1]])[[3]]
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB")))
coef <- scan1coef(prob2, phe, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("SS", "SB", "BB", "sex", "SB:sex", "BB:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(n_ind(iron)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, zerosum=FALSE)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, se=TRUE, zerosum=FALSE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")),
list(posnames, c("SS", "SB", "BB", "sex", "another", "SB:sex", "BB:sex",
"SB:another", "BB:another")))
})
test_that("scan1coef for intercross, with contrasts", {
skip_on_cran()
set.seed(9308594)
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
prob2 <- calc_genoprob(iron[,7], map, error_prob=0.002) # just look at chr 7
phe <- iron$pheno[,1] # liver phenotype
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
covar <- cbind(sex=covar)
weights <- runif(n_ind(iron), 0, 5)
names(weights) <- names(phe)
# different organization of probs for qtl2 and lm() code
prob <- aperm(prob2[[1]], c(1,3,2)) # rearrange as expected for lm()
# use contrasts
contrasts <- cbind(mu=c(1,1,1), a=c(-1, 0, 1), d=c(0, 1,0))
# change to contrasts for use with lm()
dn <- dimnames(prob)
p <- prob
dim(p) <- c(prod(dim(prob)[1:2]), dim(prob)[3])
p <- p %*% contrasts
dim(p) <- dim(prob)
prob <- p
dimnames(prob) <- list(dn[[1]], dn[[2]], colnames(contrasts))
posnames <- dimnames(prob2[[1]])[[3]]
# no covariates
expected <- eff_via_lm(prob, phe)
coef <- scan1coef(prob2, phe, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
coef <- scan1coef(prob2, phe, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d")))
# no covariates, weighted
expected <- eff_via_lm(prob, phe, weights=weights)
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
coef <- scan1coef(prob2, phe, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d")))
# one add've covariate
expected <- eff_via_lm(prob, phe, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex")))
# one add've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex")))
# one int've covariate
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
# one int've covariate, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "a:sex", "d:sex")))
# two int've covariates
covar <- cbind(covar, another=rnorm(n_ind(iron)))
expected <- eff_via_lm(prob, phe, covar, covar)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
expect_equivalent(coef, expected)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
# two int've covariates, weighted
expected <- eff_via_lm(prob, phe, covar, covar, weights=weights)
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames,
c("mu", "a", "d", "sex", "another", "a:sex", "d:sex",
"a:another", "d:another")))
coef <- scan1coef(prob2, phe, addcovar=covar, intcovar=covar, weights=weights, contrasts=contrasts, se=TRUE)
expect_equivalent(coef, expected)
expect_equal(dimnames(coef), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
expect_equivalent(attr(coef, "SE"), attr(expected, "SE"))
expect_equal(dimnames(attr(coef, "SE")), list(posnames, c("mu", "a", "d", "sex", "another",
"a:sex", "d:sex", "a:another", "d:another")))
})
test_that("scan1coef deals with mismatching individuals", {
skip_on_cran()
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
probs <- probs[,"3"]
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
phe <- iron$pheno[,2,drop=FALSE]
ind <- c(1:50, 101:150)
expected <- scan1coef(probs[ind,], phe[ind,,drop=FALSE], addcovar=X[ind])
expect_equal(scan1coef(probs[ind,], phe, addcovar=X), expected)
expect_equal(scan1coef(probs, phe[ind,,drop=FALSE], addcovar=X), expected)
expect_equal(scan1coef(probs, phe, addcovar=X[ind]), expected)
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.