Nothing
tests/testthat/test-fit1.R
In qtl2: Quantitative Trait Locus Mapping in Experimental Crosses
context("fit1")
test_that("fit1 by H-K works in intercross", {
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,c(18:19,"X")]
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, zerosum=FALSE) })
# fit1, no missing data
npos <- sapply(probs, function(a) dim(a)[3])
pmar <- c(3, 4, 12)
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, zerosum=FALSE) })
pos <- cumsum(c(0, npos[-3])) + pmar
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# repeat the whole thing with a couple of missing phenotypes
pheno[c(187, 244)] <- NA
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, se=TRUE, zerosum=FALSE) })
# fit1, missing data
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# check SEs
for(i in 1:3)
expect_equal(out_fit1[[i]]$SE, attr(coef[[i]], "SE")[pmar[i],])
# direct calculations, chr 18
lm0 <- lm(pheno ~ covar)
X <- cbind(probs[[1]][,,pmar[1]], covar)
colnames(X) <- c("SS", "SB", "BB", "ac1")
lm1 <- lm(pheno ~ -1 + X)
n <- sum(!is.na(pheno))
sigsq0 <- sum(lm0$resid^2)/n
sigsq1 <- sum(lm1$resid^2)/n
expect_equal(out_fit1[[1]]$lod, n/2*log10(sum(lm0$resid^2)/sum(lm1$resid^2)))
expect_equal(out_fit1[[1]]$ind_lod, (dnorm(lm1$resid,0,sqrt(sigsq1), TRUE) - dnorm(lm0$resid,0,sqrt(sigsq0),TRUE))/log(10))
expect_equal(out_fit1[[1]]$coef, stats::setNames(lm1$coef, c("SS", "SB", "BB", "ac1")))
expect_equal(out_fit1[[1]]$SE, stats::setNames(summary(lm1)$coef[,2], c("SS", "SB", "BB", "ac1")))
# direct calculations, chr X
lm0 <- lm(pheno ~ Xcovar)
X <- probs[[3]][,,pmar[3]]
colnames(X) <- c("SS", "SB", "BS", "BB", "SY", "BY")
lm1 <- lm(pheno ~ -1 + X)
n <- sum(!is.na(pheno))
sigsq0 <- sum(lm0$resid^2)/n
sigsq1 <- sum(lm1$resid^2)/n
expect_equal(out_fit1[[3]]$lod, n/2*log10(sum(lm0$resid^2)/sum(lm1$resid^2)))
expect_equal(out_fit1[[3]]$ind_lod, (dnorm(lm1$resid,0,sqrt(sigsq1), TRUE) - dnorm(lm0$resid,0,sqrt(sigsq0),TRUE))/log(10))
expect_equal(out_fit1[[3]]$coef, stats::setNames(lm1$coef, c("SS", "SB", "BS", "BB", "SY", "BY")))
expect_equal(out_fit1[[3]]$SE, stats::setNames(summary(lm1)$coef[,2], c("SS", "SB", "BS", "BB", "SY", "BY")))
})
test_that("fit1 by H-K works in intercross, with weights", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,c(18:19,"X")]
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
w <- setNames(runif(n_ind(iron), 1, 10), ind_ids(iron))
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar, weights=w)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, weights=w, zerosum=FALSE) })
# fit1, no missing data
npos <- sapply(probs, function(a) dim(a)[3])
pmar <- c(3, 4, 12)
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, weights=w, zerosum=FALSE) })
pos <- cumsum(c(0, npos[-3])) + pmar
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# repeat the whole thing with a couple of missing phenotypes
pheno[c(187, 244)] <- NA
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar, weights=w)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, se=TRUE, weights=w, zerosum=FALSE) })
# fit1, missing data
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# check SEs
for(i in 1:3)
expect_equal(out_fit1[[i]]$SE, attr(coef[[i]], "SE")[pmar[i],])
})
test_that("fit1 by H-K works in riself", {
grav2 <- read_cross2(system.file("extdata", "grav2.zip", package="qtl2"))
grav2 <- grav2[,4:5]
map <- insert_pseudomarkers(grav2$gmap, step=1)
probs <- calc_genoprob(grav2, map, error_prob=0.002)
pheno <- grav2$pheno[,219]
# calculate LOD scores
out <- scan1(probs, pheno)
# estimate coefficients
coef <- lapply(seq_len(length(probs)), function(i) scan1coef(subset(probs, chr=names(probs)[i]), pheno, zerosum=FALSE))
# fit1, no missing data
npos <- sapply(probs, function(a) dim(a)[3])
pmar <- c(1, 172)
out_fit1 <- lapply(seq(along=pmar), function(i) fit1(probs[[i]][,,pmar[i]], pheno, zerosum=FALSE))
pos <- c(0,npos[1]) + pmar
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:2) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:2)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# repeat the whole thing with a couple of missing phenotypes
pheno[c(24, 106)] <- NA
# calculate LOD scores
out <- scan1(probs, pheno)
# estimate coefficients
coef <- lapply(seq_len(length(probs)), function(i) scan1coef(subset(probs, chr=names(probs)[i]), pheno, se=TRUE, zerosum=FALSE))
# fit1, missing data
out_fit1 <- lapply(seq(along=pmar), function(i) fit1(probs[[i]][,,pmar[i]], pheno, se=TRUE, zerosum=FALSE))
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:2) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:2)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# check SEs
for(i in 1:2)
expect_equal(out_fit1[[i]]$SE, attr(coef[[i]], "SE")[pmar[i],])
# direct calculations, chr 18
lm0 <- lm(pheno ~ 1)
X <- probs[[1]][,,pmar[1]]
colnames(X) <- c("LL", "CC")
lm1 <- lm(pheno ~ -1 + X)
n <- sum(!is.na(pheno))
sigsq0 <- sum(lm0$resid^2)/n
sigsq1 <- sum(lm1$resid^2)/n
expect_equal(out_fit1[[1]]$lod, n/2*log10(sum(lm0$resid^2)/sum(lm1$resid^2)))
expect_equal(out_fit1[[1]]$ind_lod, (dnorm(lm1$resid,0,sqrt(sigsq1), TRUE) - dnorm(lm0$resid,0,sqrt(sigsq0),TRUE))/log(10))
expect_equal(out_fit1[[1]]$coef, stats::setNames(lm1$coef, c("LL", "CC")))
expect_equal(out_fit1[[1]]$SE, stats::setNames(summary(lm1)$coef[,2], c("LL", "CC")))
})
test_that("fit1 by LMM works in intercross", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
kinship_loco <- calc_kinship(probs, "loco")
probs <- probs[,c(7,19,"X")]
kinship_loco <- kinship_loco[c(7,19,"X")]
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
# calculate LOD scores
out <- scan1(probs, pheno, kinship, addcovar=covar, Xcovar=Xcovar)
out_loco <- scan1(probs, pheno, kinship_loco, addcovar=covar, Xcovar=Xcovar)
# estimate coefficients (with SEs)
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship,
addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
coef_loco <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship_loco[i],
addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
names(coef) <- names(coef_loco) <- names(kinship_loco)
# positions to test fit1()
chr <- rep(c(7,19,"X"), each=4)
n_pos <- sapply(probs, function(a) dim(a)[3])
index <- c(1,6,7,58, 1,15,31,36, 1,20,27,30)
pos <- index + rep(cumsum(c(0,n_pos)), c(4,4,4,0))
# fit1, no missing data
out_fit1 <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship, addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
out_fit1_loco <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship_loco[[chr[i]]], addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
# check LOD vs scan1
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1], tol=1e-6)
}
# same with _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$lod, out_loco[pos[i],1], tol=1e-6)
}
# check coefficients
for(i in seq(along=out_fit1))
expect_equal(out_fit1[[i]]$coef, coef[[chr[i]]][index[i],])
# same, _loco version
for(i in seq(along=out_fit1))
expect_equal(out_fit1_loco[[i]]$coef, coef_loco[[chr[i]]][index[i],])
# check SEs
for(i in seq(along=out_fit1))
expect_equal(out_fit1[[i]]$SE, attr(coef[[chr[i]]], "SE")[index[i],])
# same, _loco version
for(i in seq(along=out_fit1))
expect_equal(out_fit1_loco[[i]]$SE, attr(coef_loco[[chr[i]]], "SE")[index[i],])
})
test_that("fit1 by LMM works in intercross, with weights", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
kinship_loco <- calc_kinship(probs, "loco")
probs <- probs[,c(7,19,"X")]
kinship_loco <- kinship_loco[c(7,19,"X")]
set.seed(2661343)
w <- setNames(runif(n_ind(iron), 1, 8), ind_ids(iron))
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
# calculate LOD scores
out <- scan1(probs, pheno, kinship, addcovar=covar, Xcovar=Xcovar, weights=w)
out_loco <- scan1(probs, pheno, kinship_loco, addcovar=covar, Xcovar=Xcovar, weights=w)
# estimate coefficients (with SEs)
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship,
addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
coef_loco <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship_loco[i],
addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
names(coef) <- names(coef_loco) <- names(kinship_loco)
# positions to test fit1()
chr <- rep(c(7,19,"X"), each=4)
n_pos <- sapply(probs, function(a) dim(a)[3])
index <- c(1,6,7,58, 1,15,31,36, 1,20,27,30)
pos <- index + rep(cumsum(c(0,n_pos)), c(4,4,4,0))
# fit1, no missing data
out_fit1 <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship, addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
out_fit1_loco <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship_loco[[chr[i]]], addcovar=cov, nullcovar=nullcov,
se=TRUE, weights=w, zerosum=FALSE) })
# check LOD vs scan1
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1], tol=1e-6)
}
# same with _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$lod, out_loco[pos[i],1], tol=1e-6)
}
# check coefficients
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$coef, coef[[chr[i]]][index[i],])
}
# same, _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$coef, coef_loco[[chr[i]]][index[i],])
}
# check SEs
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$SE, attr(coef[[chr[i]]], "SE")[index[i],])
}
# same, _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$SE, attr(coef_loco[[chr[i]]], "SE")[index[i],])
}
})
test_that("fit1 handles contrasts properly in a backcross", {
library(qtl)
data(fake.bc)
fake.bc <- subset(fake.bc, chr=c(2,5))
fake_bc <- convert2cross2(fake.bc)
fake.bc <- calc.genoprob(fake.bc, error.prob=0.002, map.function="c-f")
pr <- calc_genoprob(fake_bc, error_prob=0.002, map_function="c-f")
# fit single-QTL model at two positions
qtl <- makeqtl(fake.bc, c(2,5), pos=c(44.8, 9.8), what="prob")
out.fq1 <- fitqtl(fake.bc, qtl=qtl, formula=y~Q1, get.ests=TRUE, method="hk")
co1 <- summary(out.fq1)$ests
out.fq2 <- fitqtl(fake.bc, qtl=qtl, formula=y~Q2, get.ests=TRUE, method="hk")
co2 <- summary(out.fq2)$ests
contrasts <- cbind(mu=c(1,1), a=c(-0.5,0.5))
out_fit1a <- fit1(pull_genoprobpos(pr, "D2M336"), fake_bc$pheno[,1],
contrasts=contrasts)
expect_equivalent(co1[,"est"], out_fit1a$coef)
expect_equivalent(co1[,"SE"], out_fit1a$SE)
out_fit1b <- fit1(pull_genoprobpos(pr, "D5M394"), fake_bc$pheno[,1],
contrasts=contrasts)
expect_equivalent(co2[,"est"], out_fit1b$coef)
expect_equivalent(co2[,"SE"], out_fit1b$SE)
# compare to scan1coef()
ests_c2 <- scan1coef(pr[,"2"], fake_bc$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c2["D2M336",], out_fit1a$coef)
expect_equal(attr(ests_c2, "SE")["D2M336",], out_fit1a$SE)
ests_c5 <- scan1coef(pr[,"5"], fake_bc$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c5["D5M394",], out_fit1b$coef)
expect_equal(attr(ests_c5, "SE")["D5M394",], out_fit1b$SE)
})
test_that("fit1 handles contrasts properly in an intercross", {
library(qtl)
data(fake.f2)
fake.f2 <- subset(fake.f2, chr=c(1,13))
fake_f2 <- convert2cross2(fake.f2)
fake.f2 <- calc.genoprob(fake.f2, error.prob=0.002, map.function="c-f")
pr <- calc_genoprob(fake_f2, error_prob=0.002, map_function="c-f")
# fit single-QTL model at two positions
qtl <- makeqtl(fake.f2, c(1,13), pos=c(37.1, 24.0), what="prob")
out.fq1 <- fitqtl(fake.f2, qtl=qtl, formula=y~Q1, get.ests=TRUE, method="hk")
co1 <- summary(out.fq1)$ests
out.fq2 <- fitqtl(fake.f2, qtl=qtl, formula=y~Q2, get.ests=TRUE, method="hk")
co2 <- summary(out.fq2)$ests
# R/qtl1 and R/qtl2 give different answers for the intercept,
# so we'll just look at the additive and dominance effects
contrasts <- cbind(mu=c(1,1,1), a=c(-1,0,1), d=c(0,1,0))
out_fit1a <- fit1(pull_genoprobpos(pr, "D1M437"), fake_f2$pheno[,1],
contrasts=contrasts)
expect_equivalent(co1[-1,"est"], out_fit1a$coef[-1])
expect_equivalent(co1[-1,"SE"], out_fit1a$SE[-1])
out_fit1b <- fit1(pull_genoprobpos(pr, "D13M254"), fake_f2$pheno[,1],
contrasts=contrasts)
expect_equivalent(co2[-1,"est"], out_fit1b$coef[-1])
expect_equivalent(co2[-1,"SE"], out_fit1b$SE[-1])
# compare to scan1coef()
ests_c1 <- scan1coef(pr[,"1"], fake_f2$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c1["D1M437",], out_fit1a$coef)
expect_equal(attr(ests_c1, "SE")["D1M437",], out_fit1a$SE)
ests_c13 <- scan1coef(pr[,"13"], fake_f2$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c13["D13M254",], out_fit1b$coef)
expect_equal(attr(ests_c13, "SE")["D13M254",], out_fit1b$SE)
})
test_that("fit1 works with blup=TRUE", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
Xcovar <- get_x_covar(iron)
kinship <- calc_kinship(probs, "loco")
# pull out a position
marker <- find_marker(map, "2", 56.8)
pr_pos <- pull_genoprobpos(probs, marker=marker)
# reduce the probs to a few positions
probs[[2]] <- probs[[2]][,, which(dimnames(probs[[2]])[[3]]==marker) + -2:1]
out_scan1blup <- scan1blup(probs[,2], iron$pheno[,1], kinship[[2]])
out_scan1blup_nok <- scan1blup(probs[,2], iron$pheno[,1])
expect_equal( fit1(pr_pos, iron$pheno[,1], kinship[[2]], blup=TRUE, se=FALSE)$coef,
out_scan1blup[marker,] )
expect_equal( fit1(pr_pos, iron$pheno[,1], blup=TRUE, se=FALSE)$coef,
out_scan1blup_nok[marker,] )
})
test_that("fit1 fitted values don't depend on order of individuals", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
out <- scan1(probs[,"7"], pheno, addcovar=covar)
max_pos <- max(out, map)
pr_max <- pull_genoprobpos(probs, map, max_pos$chr, max_pos$pos)
out_fit1 <- fit1(pr_max, pheno, addcovar=covar)
ind <- rownames(pr_max)
pr_max_samp <- pr_max[sample(ind),]
out_fit1b <- fit1(pr_max_samp, pheno, addcovar=covar)
testthat::expect_equal(out_fit1$fitted[ind] , out_fit1b$fitted[ind])
k <- calc_kinship(probs)
out_fit1 <- fit1(pr_max, pheno, addcovar=covar, kinship=k)
out_fit1b <- fit1(pr_max_samp[ind,], pheno[ind], addcovar=covar[ind], kinship=k[ind,ind])
testthat::expect_equal(out_fit1$fitted[ind] , out_fit1b$fitted[ind])
})
test_that("fit1 works without genoprobs", {
set.seed(20201215)
n <- 100
nam <- paste0("mouse", sample(10*n, n))
phe <- setNames(rnorm(n), nam)
cov <- setNames(sample(0:1, n, replace=TRUE), nam)
lm_out <- lm(phe ~ cov)
lm_sum <- summary(lm_out)
coef_names <- c("intercept", "ac1", "intercept")
expected <- list(lod=0,
ind_lod=setNames(rep(0, n), nam),
coef=setNames(c(0, lm_out$coef[2], lm_out$coef[1]), coef_names),
SE=setNames(lm_sum$coef[c(1,2,1),"Std. Error"], coef_names),
fitted=lm_out$fitted,
resid=lm_out$resid)
expect_equal(fit1(pheno=phe, addcovar=cov), expected)
k <- matrix(0.5,ncol=n, nrow=n)
diag(k) <- 1
dimnames(k) <- list(nam, nam)
# just test that this works
should_work <- fit1(pheno=phe, addcovar=cov, kinship=k)
})
Try the qtl2 package in your browser
Any scripts or data that you put into this service are public.
qtl2 documentation built on April 22, 2023, 1:10 a.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("fit1")
test_that("fit1 by H-K works in intercross", {
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,c(18:19,"X")]
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, zerosum=FALSE) })
# fit1, no missing data
npos <- sapply(probs, function(a) dim(a)[3])
pmar <- c(3, 4, 12)
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, zerosum=FALSE) })
pos <- cumsum(c(0, npos[-3])) + pmar
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# repeat the whole thing with a couple of missing phenotypes
pheno[c(187, 244)] <- NA
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, se=TRUE, zerosum=FALSE) })
# fit1, missing data
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# check SEs
for(i in 1:3)
expect_equal(out_fit1[[i]]$SE, attr(coef[[i]], "SE")[pmar[i],])
# direct calculations, chr 18
lm0 <- lm(pheno ~ covar)
X <- cbind(probs[[1]][,,pmar[1]], covar)
colnames(X) <- c("SS", "SB", "BB", "ac1")
lm1 <- lm(pheno ~ -1 + X)
n <- sum(!is.na(pheno))
sigsq0 <- sum(lm0$resid^2)/n
sigsq1 <- sum(lm1$resid^2)/n
expect_equal(out_fit1[[1]]$lod, n/2*log10(sum(lm0$resid^2)/sum(lm1$resid^2)))
expect_equal(out_fit1[[1]]$ind_lod, (dnorm(lm1$resid,0,sqrt(sigsq1), TRUE) - dnorm(lm0$resid,0,sqrt(sigsq0),TRUE))/log(10))
expect_equal(out_fit1[[1]]$coef, stats::setNames(lm1$coef, c("SS", "SB", "BB", "ac1")))
expect_equal(out_fit1[[1]]$SE, stats::setNames(summary(lm1)$coef[,2], c("SS", "SB", "BB", "ac1")))
# direct calculations, chr X
lm0 <- lm(pheno ~ Xcovar)
X <- probs[[3]][,,pmar[3]]
colnames(X) <- c("SS", "SB", "BS", "BB", "SY", "BY")
lm1 <- lm(pheno ~ -1 + X)
n <- sum(!is.na(pheno))
sigsq0 <- sum(lm0$resid^2)/n
sigsq1 <- sum(lm1$resid^2)/n
expect_equal(out_fit1[[3]]$lod, n/2*log10(sum(lm0$resid^2)/sum(lm1$resid^2)))
expect_equal(out_fit1[[3]]$ind_lod, (dnorm(lm1$resid,0,sqrt(sigsq1), TRUE) - dnorm(lm0$resid,0,sqrt(sigsq0),TRUE))/log(10))
expect_equal(out_fit1[[3]]$coef, stats::setNames(lm1$coef, c("SS", "SB", "BS", "BB", "SY", "BY")))
expect_equal(out_fit1[[3]]$SE, stats::setNames(summary(lm1)$coef[,2], c("SS", "SB", "BS", "BB", "SY", "BY")))
})
test_that("fit1 by H-K works in intercross, with weights", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,c(18:19,"X")]
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
w <- setNames(runif(n_ind(iron), 1, 10), ind_ids(iron))
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar, weights=w)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, weights=w, zerosum=FALSE) })
# fit1, no missing data
npos <- sapply(probs, function(a) dim(a)[3])
pmar <- c(3, 4, 12)
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, weights=w, zerosum=FALSE) })
pos <- cumsum(c(0, npos[-3])) + pmar
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# repeat the whole thing with a couple of missing phenotypes
pheno[c(187, 244)] <- NA
# calculate LOD scores
out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar, weights=w)
# estimate coefficients; no covariates for X chromosome
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) cov <- NULL
else cov <- covar
scan1coef(subset(probs, chr=names(probs)[i]), pheno, addcovar=cov, se=TRUE, weights=w, zerosum=FALSE) })
# fit1, missing data
out_fit1 <- lapply(seq(along=pmar),
function(i) {
if(i==3) { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[i]][,,pmar[i]], pheno, addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:3) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:3)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# check SEs
for(i in 1:3)
expect_equal(out_fit1[[i]]$SE, attr(coef[[i]], "SE")[pmar[i],])
})
test_that("fit1 by H-K works in riself", {
grav2 <- read_cross2(system.file("extdata", "grav2.zip", package="qtl2"))
grav2 <- grav2[,4:5]
map <- insert_pseudomarkers(grav2$gmap, step=1)
probs <- calc_genoprob(grav2, map, error_prob=0.002)
pheno <- grav2$pheno[,219]
# calculate LOD scores
out <- scan1(probs, pheno)
# estimate coefficients
coef <- lapply(seq_len(length(probs)), function(i) scan1coef(subset(probs, chr=names(probs)[i]), pheno, zerosum=FALSE))
# fit1, no missing data
npos <- sapply(probs, function(a) dim(a)[3])
pmar <- c(1, 172)
out_fit1 <- lapply(seq(along=pmar), function(i) fit1(probs[[i]][,,pmar[i]], pheno, zerosum=FALSE))
pos <- c(0,npos[1]) + pmar
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:2) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:2)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# repeat the whole thing with a couple of missing phenotypes
pheno[c(24, 106)] <- NA
# calculate LOD scores
out <- scan1(probs, pheno)
# estimate coefficients
coef <- lapply(seq_len(length(probs)), function(i) scan1coef(subset(probs, chr=names(probs)[i]), pheno, se=TRUE, zerosum=FALSE))
# fit1, missing data
out_fit1 <- lapply(seq(along=pmar), function(i) fit1(probs[[i]][,,pmar[i]], pheno, se=TRUE, zerosum=FALSE))
# check LOD vs scan1, plus ind'l contributions to LOD
for(i in 1:2) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1])
expect_equal(sum(out_fit1[[i]]$ind_lod), out_fit1[[i]]$lod)
}
# check coefficients
for(i in 1:2)
expect_equal(out_fit1[[i]]$coef, coef[[i]][pmar[i],])
# check SEs
for(i in 1:2)
expect_equal(out_fit1[[i]]$SE, attr(coef[[i]], "SE")[pmar[i],])
# direct calculations, chr 18
lm0 <- lm(pheno ~ 1)
X <- probs[[1]][,,pmar[1]]
colnames(X) <- c("LL", "CC")
lm1 <- lm(pheno ~ -1 + X)
n <- sum(!is.na(pheno))
sigsq0 <- sum(lm0$resid^2)/n
sigsq1 <- sum(lm1$resid^2)/n
expect_equal(out_fit1[[1]]$lod, n/2*log10(sum(lm0$resid^2)/sum(lm1$resid^2)))
expect_equal(out_fit1[[1]]$ind_lod, (dnorm(lm1$resid,0,sqrt(sigsq1), TRUE) - dnorm(lm0$resid,0,sqrt(sigsq0),TRUE))/log(10))
expect_equal(out_fit1[[1]]$coef, stats::setNames(lm1$coef, c("LL", "CC")))
expect_equal(out_fit1[[1]]$SE, stats::setNames(summary(lm1)$coef[,2], c("LL", "CC")))
})
test_that("fit1 by LMM works in intercross", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
kinship_loco <- calc_kinship(probs, "loco")
probs <- probs[,c(7,19,"X")]
kinship_loco <- kinship_loco[c(7,19,"X")]
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
# calculate LOD scores
out <- scan1(probs, pheno, kinship, addcovar=covar, Xcovar=Xcovar)
out_loco <- scan1(probs, pheno, kinship_loco, addcovar=covar, Xcovar=Xcovar)
# estimate coefficients (with SEs)
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship,
addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
coef_loco <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship_loco[i],
addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
names(coef) <- names(coef_loco) <- names(kinship_loco)
# positions to test fit1()
chr <- rep(c(7,19,"X"), each=4)
n_pos <- sapply(probs, function(a) dim(a)[3])
index <- c(1,6,7,58, 1,15,31,36, 1,20,27,30)
pos <- index + rep(cumsum(c(0,n_pos)), c(4,4,4,0))
# fit1, no missing data
out_fit1 <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship, addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
out_fit1_loco <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship_loco[[chr[i]]], addcovar=cov, nullcovar=nullcov, se=TRUE, zerosum=FALSE) })
# check LOD vs scan1
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1], tol=1e-6)
}
# same with _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$lod, out_loco[pos[i],1], tol=1e-6)
}
# check coefficients
for(i in seq(along=out_fit1))
expect_equal(out_fit1[[i]]$coef, coef[[chr[i]]][index[i],])
# same, _loco version
for(i in seq(along=out_fit1))
expect_equal(out_fit1_loco[[i]]$coef, coef_loco[[chr[i]]][index[i],])
# check SEs
for(i in seq(along=out_fit1))
expect_equal(out_fit1[[i]]$SE, attr(coef[[chr[i]]], "SE")[index[i],])
# same, _loco version
for(i in seq(along=out_fit1))
expect_equal(out_fit1_loco[[i]]$SE, attr(coef_loco[[chr[i]]], "SE")[index[i],])
})
test_that("fit1 by LMM works in intercross, with weights", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
kinship_loco <- calc_kinship(probs, "loco")
probs <- probs[,c(7,19,"X")]
kinship_loco <- kinship_loco[c(7,19,"X")]
set.seed(2661343)
w <- setNames(runif(n_ind(iron), 1, 8), ind_ids(iron))
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
Xcovar <- get_x_covar(iron)
# calculate LOD scores
out <- scan1(probs, pheno, kinship, addcovar=covar, Xcovar=Xcovar, weights=w)
out_loco <- scan1(probs, pheno, kinship_loco, addcovar=covar, Xcovar=Xcovar, weights=w)
# estimate coefficients (with SEs)
coef <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship,
addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
coef_loco <- lapply(seq_len(length(probs)), function(i) {
if(i==3) { cov <- NULL; nullcov <- Xcovar }
else { cov <- covar; nullcov <- NULL }
scan1coef(subset(probs, chr=names(probs)[i]), pheno, kinship_loco[i],
addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
names(coef) <- names(coef_loco) <- names(kinship_loco)
# positions to test fit1()
chr <- rep(c(7,19,"X"), each=4)
n_pos <- sapply(probs, function(a) dim(a)[3])
index <- c(1,6,7,58, 1,15,31,36, 1,20,27,30)
pos <- index + rep(cumsum(c(0,n_pos)), c(4,4,4,0))
# fit1, no missing data
out_fit1 <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship, addcovar=cov, nullcovar=nullcov, se=TRUE, weights=w, zerosum=FALSE) })
out_fit1_loco <- lapply(seq(along=chr),
function(i) {
if(chr[i]=="X") { nullcov <- Xcovar; cov <- NULL } # need Xcovar under null on X chr but no other covariates
else { nullcov <- NULL; cov <- covar } # sex as covariate; no additional covariates under null
fit1(probs[[chr[i]]][,,index[i]], pheno, kinship_loco[[chr[i]]], addcovar=cov, nullcovar=nullcov,
se=TRUE, weights=w, zerosum=FALSE) })
# check LOD vs scan1
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$lod, out[pos[i],1], tol=1e-6)
}
# same with _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$lod, out_loco[pos[i],1], tol=1e-6)
}
# check coefficients
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$coef, coef[[chr[i]]][index[i],])
}
# same, _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$coef, coef_loco[[chr[i]]][index[i],])
}
# check SEs
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1[[i]]$SE, attr(coef[[chr[i]]], "SE")[index[i],])
}
# same, _loco version
for(i in seq(along=out_fit1)) {
expect_equal(out_fit1_loco[[i]]$SE, attr(coef_loco[[chr[i]]], "SE")[index[i],])
}
})
test_that("fit1 handles contrasts properly in a backcross", {
library(qtl)
data(fake.bc)
fake.bc <- subset(fake.bc, chr=c(2,5))
fake_bc <- convert2cross2(fake.bc)
fake.bc <- calc.genoprob(fake.bc, error.prob=0.002, map.function="c-f")
pr <- calc_genoprob(fake_bc, error_prob=0.002, map_function="c-f")
# fit single-QTL model at two positions
qtl <- makeqtl(fake.bc, c(2,5), pos=c(44.8, 9.8), what="prob")
out.fq1 <- fitqtl(fake.bc, qtl=qtl, formula=y~Q1, get.ests=TRUE, method="hk")
co1 <- summary(out.fq1)$ests
out.fq2 <- fitqtl(fake.bc, qtl=qtl, formula=y~Q2, get.ests=TRUE, method="hk")
co2 <- summary(out.fq2)$ests
contrasts <- cbind(mu=c(1,1), a=c(-0.5,0.5))
out_fit1a <- fit1(pull_genoprobpos(pr, "D2M336"), fake_bc$pheno[,1],
contrasts=contrasts)
expect_equivalent(co1[,"est"], out_fit1a$coef)
expect_equivalent(co1[,"SE"], out_fit1a$SE)
out_fit1b <- fit1(pull_genoprobpos(pr, "D5M394"), fake_bc$pheno[,1],
contrasts=contrasts)
expect_equivalent(co2[,"est"], out_fit1b$coef)
expect_equivalent(co2[,"SE"], out_fit1b$SE)
# compare to scan1coef()
ests_c2 <- scan1coef(pr[,"2"], fake_bc$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c2["D2M336",], out_fit1a$coef)
expect_equal(attr(ests_c2, "SE")["D2M336",], out_fit1a$SE)
ests_c5 <- scan1coef(pr[,"5"], fake_bc$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c5["D5M394",], out_fit1b$coef)
expect_equal(attr(ests_c5, "SE")["D5M394",], out_fit1b$SE)
})
test_that("fit1 handles contrasts properly in an intercross", {
library(qtl)
data(fake.f2)
fake.f2 <- subset(fake.f2, chr=c(1,13))
fake_f2 <- convert2cross2(fake.f2)
fake.f2 <- calc.genoprob(fake.f2, error.prob=0.002, map.function="c-f")
pr <- calc_genoprob(fake_f2, error_prob=0.002, map_function="c-f")
# fit single-QTL model at two positions
qtl <- makeqtl(fake.f2, c(1,13), pos=c(37.1, 24.0), what="prob")
out.fq1 <- fitqtl(fake.f2, qtl=qtl, formula=y~Q1, get.ests=TRUE, method="hk")
co1 <- summary(out.fq1)$ests
out.fq2 <- fitqtl(fake.f2, qtl=qtl, formula=y~Q2, get.ests=TRUE, method="hk")
co2 <- summary(out.fq2)$ests
# R/qtl1 and R/qtl2 give different answers for the intercept,
# so we'll just look at the additive and dominance effects
contrasts <- cbind(mu=c(1,1,1), a=c(-1,0,1), d=c(0,1,0))
out_fit1a <- fit1(pull_genoprobpos(pr, "D1M437"), fake_f2$pheno[,1],
contrasts=contrasts)
expect_equivalent(co1[-1,"est"], out_fit1a$coef[-1])
expect_equivalent(co1[-1,"SE"], out_fit1a$SE[-1])
out_fit1b <- fit1(pull_genoprobpos(pr, "D13M254"), fake_f2$pheno[,1],
contrasts=contrasts)
expect_equivalent(co2[-1,"est"], out_fit1b$coef[-1])
expect_equivalent(co2[-1,"SE"], out_fit1b$SE[-1])
# compare to scan1coef()
ests_c1 <- scan1coef(pr[,"1"], fake_f2$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c1["D1M437",], out_fit1a$coef)
expect_equal(attr(ests_c1, "SE")["D1M437",], out_fit1a$SE)
ests_c13 <- scan1coef(pr[,"13"], fake_f2$pheno[,1], contrasts=contrasts, se=TRUE)
expect_equal(ests_c13["D13M254",], out_fit1b$coef)
expect_equal(attr(ests_c13, "SE")["D13M254",], out_fit1b$SE)
})
test_that("fit1 works with blup=TRUE", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=1)
probs <- calc_genoprob(iron, map, error_prob=0.002)
Xcovar <- get_x_covar(iron)
kinship <- calc_kinship(probs, "loco")
# pull out a position
marker <- find_marker(map, "2", 56.8)
pr_pos <- pull_genoprobpos(probs, marker=marker)
# reduce the probs to a few positions
probs[[2]] <- probs[[2]][,, which(dimnames(probs[[2]])[[3]]==marker) + -2:1]
out_scan1blup <- scan1blup(probs[,2], iron$pheno[,1], kinship[[2]])
out_scan1blup_nok <- scan1blup(probs[,2], iron$pheno[,1])
expect_equal( fit1(pr_pos, iron$pheno[,1], kinship[[2]], blup=TRUE, se=FALSE)$coef,
out_scan1blup[marker,] )
expect_equal( fit1(pr_pos, iron$pheno[,1], blup=TRUE, se=FALSE)$coef,
out_scan1blup_nok[marker,] )
})
test_that("fit1 fitted values don't depend on order of individuals", {
skip_if(isnt_karl(), "this test only run locally")
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
map <- insert_pseudomarkers(iron$gmap, step=5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
pheno <- iron$pheno[,1]
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)
out <- scan1(probs[,"7"], pheno, addcovar=covar)
max_pos <- max(out, map)
pr_max <- pull_genoprobpos(probs, map, max_pos$chr, max_pos$pos)
out_fit1 <- fit1(pr_max, pheno, addcovar=covar)
ind <- rownames(pr_max)
pr_max_samp <- pr_max[sample(ind),]
out_fit1b <- fit1(pr_max_samp, pheno, addcovar=covar)
testthat::expect_equal(out_fit1$fitted[ind] , out_fit1b$fitted[ind])
k <- calc_kinship(probs)
out_fit1 <- fit1(pr_max, pheno, addcovar=covar, kinship=k)
out_fit1b <- fit1(pr_max_samp[ind,], pheno[ind], addcovar=covar[ind], kinship=k[ind,ind])
testthat::expect_equal(out_fit1$fitted[ind] , out_fit1b$fitted[ind])
})
test_that("fit1 works without genoprobs", {
set.seed(20201215)
n <- 100
nam <- paste0("mouse", sample(10*n, n))
phe <- setNames(rnorm(n), nam)
cov <- setNames(sample(0:1, n, replace=TRUE), nam)
lm_out <- lm(phe ~ cov)
lm_sum <- summary(lm_out)
coef_names <- c("intercept", "ac1", "intercept")
expected <- list(lod=0,
ind_lod=setNames(rep(0, n), nam),
coef=setNames(c(0, lm_out$coef[2], lm_out$coef[1]), coef_names),
SE=setNames(lm_sum$coef[c(1,2,1),"Std. Error"], coef_names),
fitted=lm_out$fitted,
resid=lm_out$resid)
expect_equal(fit1(pheno=phe, addcovar=cov), expected)
k <- matrix(0.5,ncol=n, nrow=n)
diag(k) <- 1
dimnames(k) <- list(nam, nam)
# just test that this works
should_work <- fit1(pheno=phe, addcovar=cov, kinship=k)
})
Try the qtl2 package in your browser
Any scripts or data that you put into this service are public.
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.