tests/testthat/test-scan1_pg.R
In rqtl/qtl2: Quantitative Trait Locus Mapping in Experimental Crosses
context("LMM genome scan by scan1 with kinship matrix")
test_that("scan1 with kinship with intercross, vs ported lmmlite code", {
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[1:10,c(1,3)]
map <- insert_pseudomarkers(iron$gmap, step=5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
out_reml <- scan1(probs, iron$pheno, kinship, reml=TRUE)
out_ml <- scan1(probs, iron$pheno, kinship, reml=FALSE)
# "by hand" calculation
y <- iron$pheno
X <- cbind(rep(1, nrow(iron$pheno)))
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% X
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
byhand1_reml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[1,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tol=1e-7)
expect_equal(as.numeric(attr(out_ml, "hsq")[1,]),
c(byhand1_ml$hsq, byhand2_ml$hsq), tol=1e-7)
# compare chromosome 1 LOD scores
d <- dim(probs[[1]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(X, probs[[1]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[1:d,1], lod_reml1, tol=1e-7)
expect_equal(out_reml[1:d,2], lod_reml2, tol=1e-7)
expect_equal(out_ml[1:d,1], lod_ml1, tol=1e-7)
expect_equal(out_ml[1:d,2], lod_ml2, tol=1e-7)
})
test_that("scan1 with intercross with X covariates for null", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
Xc <- get_x_covar(iron)
out_reml <- scan1(probs, iron$pheno, kinship, Xcovar=Xc, reml=TRUE)
out_ml <- scan1(probs, iron$pheno, kinship, Xcovar=Xc, reml=FALSE)
# "by hand" calculation
y <- iron$pheno
X <- cbind(rep(1, nrow(iron$pheno)))
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% X
Xcp <- Ke$vectors %*% Xc
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
byhand1_reml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp), reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp), reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp), reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp), reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[2,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tolerance=1e-6)
expect_equal(as.numeric(attr(out_ml, "hsq")[2,]),
c(byhand1_ml$hsq, byhand2_ml$hsq))
# compare chromosome X LOD scores
d <- dim(probs[["X"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, probs[["X"]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
index <- nrow(out_reml) - rev(1:d) + 1
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1, tolerance=1e-6)
expect_equal(out_reml[index,2], lod_reml2, tolerance=1e-6)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
})
test_that("scan1 with kinship with intercross with an additive covariate", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc, reml=TRUE, tol=1e-12)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc, reml=FALSE, tol=1e-12)
# "by hand" calculation
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% cbind(1, X)
Xcp <- Ke$vectors %*% Xc
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
# autosome null
byhand1A_reml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2A_reml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1A_ml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2A_ml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[1,]),
c(byhand1A_reml$hsq, byhand2A_reml$hsq))
expect_equal(as.numeric(attr(out_ml, "hsq")[1,]),
c(byhand1A_ml$hsq, byhand2A_ml$hsq))
# X chr null
byhand1X_reml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand2X_reml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand1X_ml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
byhand2X_ml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[2,]),
c(byhand1X_reml$hsq, byhand2X_reml$hsq), tolerance=1e-6)
expect_equal(as.numeric(attr(out_ml, "hsq")[2,]),
c(byhand1X_ml$hsq, byhand2X_ml$hsq))
# compare chromosome 2 LOD scores
d <- dim(probs[["2"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, X, probs[["2"]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1A_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2A_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1A_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2A_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1A_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2A_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1A_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2A_ml$loglik)/log(10)
index <- dim(probs[["1"]])[3] + 1:dim(probs[["2"]])[3]
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1)
expect_equal(out_reml[index,2], lod_reml2)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
# compare chromosome X LOD scores
d <- dim(probs[["X"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, probs[["X"]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1X_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2X_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1X_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2X_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1X_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2X_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1X_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2X_ml$loglik)/log(10)
index <- nrow(out_reml) - rev(1:d) + 1
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
## FIX_ME
## REML not yet working on X chromosome, when (X, probs) is not full rank
# expect_equal(out_reml[index,1], lod_reml1)
# expect_equal(out_reml[index,2], lod_reml2)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
})
test_that("scan1 with kinship with intercross with an interactive covariate", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12, intcovar_method="lowmem")
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12, intcovar_method="lowmem")
out_reml_himem <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12, intcovar_method="highmem")
out_ml_himem <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12, intcovar_method="highmem")
# same result using "highmem" and "lowmem" methods
expect_equal(out_reml_himem, out_reml)
expect_equal(out_ml_himem, out_ml)
# "by hand" calculation
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% cbind(1, X)
Xcp <- Ke$vectors %*% Xc
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
# autosome null (same as w/o interactive covariate)
byhand1A_reml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2A_reml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1A_ml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2A_ml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[1,]),
c(byhand1A_reml$hsq, byhand2A_reml$hsq))
expect_equal(as.numeric(attr(out_ml, "hsq")[1,]),
c(byhand1A_ml$hsq, byhand2A_ml$hsq))
# X chr null (same as w/o interactive covariate)
byhand1X_reml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand2X_reml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand1X_ml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
byhand2X_ml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[2,]),
c(byhand1X_reml$hsq, byhand2X_reml$hsq), tolerance=1e-6)
expect_equal(as.numeric(attr(out_ml, "hsq")[2,]),
c(byhand1X_ml$hsq, byhand2X_ml$hsq))
# compare chromosome 4 LOD scores
npos <- sapply(probs, function(a) dim(a)[[3]])
d <- npos["4"]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, X, probs[["4"]][,-1,i], probs[["4"]][,-1,i]*X)
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1A_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2A_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1A_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2A_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1A_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2A_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1A_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2A_ml$loglik)/log(10)
index <- sum(npos[1:3]) + 1:npos[4]
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1)
expect_equal(out_reml[index,2], lod_reml2)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
# compare chromosome X LOD scores
d <- dim(probs[["X"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, 3)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, X, probs[["X"]][,-1,i], probs[["X"]][,-1,i]*X)
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1X_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2X_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1X_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2X_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1X_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2X_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1X_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2X_ml$loglik)/log(10)
index <- nrow(out_reml) - rev(1:d) + 1
## FIX ME
## Not yet working on X chromosome, when (X, probs) is not full rank
# expect_equal(out_reml[index,1], lod_reml1)
# expect_equal(out_reml[index,2], lod_reml2)
# expect_equal(out_ml[index,1], lod_ml1)
# expect_equal(out_ml[index,2], lod_ml2)
})
test_that("scan1 with kinship works with LOCO, additive covariates", {
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)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
# double the eigenvalues (== kinship matrix * 2)
Ke <- lapply(Ke, function(a) { a$values <- 2*a$values; a})
# compare chromosomes 1, 6, 9, 18
chrs <- paste(c(1,6,9,18))
npos <- sapply(probs, function(a) dim(a)[[3]])
for(chr in chrs) {
nchr <- which(names(npos) == chr)
d <- npos[chr]
yp <- Ke[[chr]]$vectors %*% y
Xp <- Ke[[chr]]$vectors %*% cbind(1, X)
# autosome null
byhand1_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[nchr,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tolerance=1e-5)
expect_equal(as.numeric(attr(out_ml, "hsq")[nchr,]),
c(byhand1_ml$hsq, byhand2_ml$hsq), tolerance=1e-6)
# chromosome scan
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke[[chr]]$vectors %*% cbind(1, X, probs[[chr]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
if(nchr > 1) index <- sum(npos[1:(nchr-1)]) + 1:d
else index <- 1:d
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1, tolerance=1e-6)
expect_equal(out_reml[index,2], lod_reml2, tolerance=1e-5)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
}
})
test_that("scan1 with kinship works with LOCO, interactive covariates", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
# double the eigenvalues (== kinship matrix * 2)
Ke <- lapply(Ke, function(a) { a$values <- 2*a$values; a})
# compare chromosomes 1, 6, 9, 18
chrs <- paste(c(1,6,9,18))
npos <- sapply(probs, function(a) dim(a)[[3]])
for(chr in chrs) {
nchr <- which(names(npos) == chr)
d <- npos[chr]
yp <- Ke[[chr]]$vectors %*% y
Xp <- Ke[[chr]]$vectors %*% cbind(1, X)
# autosome null (same as w/o interactive covariate)
byhand1_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[nchr,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tolerance=1e-5)
expect_equal(as.numeric(attr(out_ml, "hsq")[nchr,]),
c(byhand1_ml$hsq, byhand2_ml$hsq), tolerance=1e-6)
# chromosome scan
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke[[chr]]$vectors %*% cbind(1, X, probs[[chr]][,-1,i], probs[[chr]][,-1,i]*X)
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
if(nchr > 1) index <- sum(npos[1:(nchr-1)]) + 1:d
else index <- 1:d
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1)
expect_equal(out_reml[index,2], lod_reml2, tolerance=1e-6)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
}
})
test_that("scan1 with kinship works with multicore", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_reml_4core <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12, cores=2)
expect_equal(out_reml, out_reml_4core)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
out_ml_4core <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12, cores=2)
expect_equal(out_ml, out_ml_4core)
})
test_that("scan1 with kinship LOD results invariant to change in scale to pheno and covar", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_reml_scale <- scan1(probs, iron$pheno/100, kinship, addcovar=X*2, intcovar=X*2,
Xcovar=Xc*2, reml=TRUE, tol=1e-12)
expect_equal(out_reml, out_reml_scale, tol=1e-6)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
out_ml_scale <- scan1(probs, iron$pheno/100, kinship, addcovar=X*4, intcovar=X*4,
Xcovar=Xc*4, reml=FALSE, tol=1e-12)
expect_equal(out_ml, out_ml_scale, tol=1e-6)
})
test_that("scan1 deals with mismatching 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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
ind <- c(1:50, 101:150)
subK <- lapply(kinship, "[", ind, ind)
expected <- scan1(probs[ind,], iron$pheno[ind,,drop=FALSE], subK, addcovar=X[ind], intcovar=X[ind],
Xcovar=Xc[ind,], reml=TRUE, tol=1e-12)
expect_equal(scan1(probs[ind,], iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno[ind,], kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, subK, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X[ind], intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X[ind],
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc[ind,], reml=TRUE, tol=1e-12), expected)
})
test_that("scan1 with weights and kinship", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
chr <- c("7", "12", "X")
probs <- probs[,chr]
kinship <- kinship[chr]
RNGkind("Mersenne-Twister") # make sure we're using the standard RNG
set.seed(28915967)
weights <- stats::setNames(sample(1:10, n_ind(iron), replace=TRUE), ind_ids(iron))
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X,
Xcovar=Xc, weights=weights, reml=TRUE, tol=1e-12)
### results via regress library
# library(regress)
# hsq <- matrix(nrow=3, ncol=2)
# dimnames(hsq) <- list(c("7", "12", "X"), c("liver", "spleen"))
# for(i in 1:2) {
# for(j in 1:3) {
# k <- kinship[[j]]*2
# if(j==3) co <- Xc else co <- X
# out_regress <- regress(iron$pheno[,i] ~ co, ~ k + diag(1/weights), identity=FALSE, tol=1e-8)
# sig <- out_regress$sigma
# hsq[j,i] <- sig[1]/sum(sig)
# }
# }
hsq <- structure(c(0.281878438466624, 0.353558919335662, 0.354315217715963,
0.20945078570429, 0.217671918799751, 0.201682575651958),
.Dim = 3:2, .Dimnames = list(c("7", "12", "X"), c("liver", "spleen")))
expect_equal(attr(out_reml, "hsq"), hsq, tol=1e-5)
### calculate lod scors from scan
# hsq <- attr(out_reml, "hsq")
# result <- vector("list", 6)
# for(i in 1:2) {
# for(j in 1:3) {
# # cholesky decomp of kinhip matrix
# v <- hsq[j,i]*2*kinship[[j]] + (1-hsq[j,i])*diag(1/weights)
# d <- solve(t(chol( v )))
# # pre-multiply phenotype and covar by t(d)
# y <- d %*% iron$pheno[,i]
# if(j==3) co <- Xc else co <- X
# co <- d %*% cbind(1, co)
#
# # log lik under null
# r <- lm(y ~ -1 + co)$resid
# llik0 <- sum(dnorm(r, 0, sqrt(mean(r^2)), log=TRUE))
#
# llik1 <- apply(probs[[j]], 3, function(z) {
# thisX <- cbind(co, d %*% z[,-1,drop=FALSE])
# r <- lm(y ~ -1 + thisX)$resid
# sum(dnorm(r, 0, sqrt(mean(r^2)), log=TRUE)) })
#
# result[[(i-1)*3+j]] <- (llik1 - llik0)/log(10)
#
# }
# }
# result <- rbind(cbind(result[[1]], result[[4]]),
# cbind(result[[2]], result[[5]]),
# cbind(result[[3]], result[[6]]))
# dimnames(result) <- dimnames(out_reml)
# class(result) <- class(out_reml)
result <- structure(c(0.765379546370331, 1.33697830149196, 2.05141217249269,
2.82748389184348, 3.57825928288821, 4.11482576884623,
4.18208695088672, 4.5010472223741, 4.77640847444171,
4.98121670508953, 5.07224467963965, 5.28435499128555,
5.10772158608095, 3.83479814287182, 3.89931966500793,
4.24832760135268, 4.18515409936768, 4.61791967696627,
4.70308738505572, 4.58758777783123, 4.88425659502959,
5.39302146958257, 5.86125538183219, 6.27973204425698,
6.62245498126028, 6.84349072896544, 6.91418168584823,
0.492648140597935, 0.483021280268028, 0.465046322957041,
0.435712712087884, 0.391979816567206, 0.332013881699469,
0.257330299826467, 0.174985513110401, 0.0975448556296765,
0.0390292238934769, 0.0087954970835623, 0.0081071126160765,
0.0315063370227757, 0.0706853835370534, 0.117761588763882,
0.166806834360401, 0.180329873353274, 0.146302909441694,
0.161137797908763, 0.192675853649471, 0.243358929100302,
0.313640413846475, 0.401241925760445, 0.501210880937258,
0.606957139785506, 0.711817493234455, 0.810381904299065,
0.899112752704855, 0.976303366208711, 1.00118015673376,
0.68879116519848, 1.13593060562675, 1.59214070065802,
1.93410070637396, 2.08715532129677, 2.08161469298463,
2.11418090796465, 2.24129850312586, 2.28750006399176,
2.24152015532401, 2.15435105671462, 2.16059659637594,
1.40659065214545, 0.389546568631669, 0.392567989473253,
0.432129729053773, 0.445384879284572, 0.468919719042007,
0.561846602500791, 0.61487493505232, 0.680306211651546,
0.828892839652596, 1.03747954481922, 1.3178288863555,
1.64752359905175, 1.97570834069646, 2.25872692788406,
1.76552836513385, 1.9624748694232, 2.179534029569,
2.40709075366897, 2.62706355895772, 2.81242765604838,
2.93166656075001, 2.95963946310072, 2.89030386259124,
2.74110806468527, 2.54394356128059, 2.33040574729971,
2.12258054254516, 1.93206224497564, 1.76308130810039,
1.61583945469733, 1.57830829020771, 1.91298146944818,
1.89937590028851, 1.86320092206996, 1.79734091344691,
1.69690361501555, 1.56240479342725, 1.4020247631517,
1.23069582035597, 1.06561266027339, 0.920732232061234,
0.803561188099648, 0.715279174632663, 0.690066084958058),
.Dim = c(57L, 2L),
.Dimnames = list(c("D7Mit74", "c7.loc4", "c7.loc6", "c7.loc9",
"c7.loc11", "D7Mit25", "c7.loc14", "c7.loc16",
"c7.loc19", "c7.loc21", "D7Nds5", "c7.loc24",
"c7.loc26", "D7mit30", "c7.loc29", "c7.loc31",
"D7Mit31", "c7.loc34", "c7.loc36", "D7Mit17",
"c7.loc39", "c7.loc41", "c7.loc44", "c7.loc46",
"c7.loc49", "c7.loc51", "D7Mit71", "D12Mit88",
"c12.loc22", "c12.loc25", "c12.loc27", "c12.loc30",
"c12.loc32", "c12.loc35", "c12.loc37", "c12.loc40",
"c12.loc42", "c12.loc45", "c12.loc47", "c12.loc50",
"c12.loc52", "c12.loc55", "c12.loc57", "D12Mit134",
"DXMit16", "cX.loc32", "cX.loc34", "cX.loc37",
"cX.loc40", "cX.loc42", "cX.loc44", "cX.loc47",
"cX.loc50", "cX.loc52", "cX.loc54",
"cX.loc57", "DXMit186"), c("liver", "spleen")),
class = c("scan1", "matrix"))
for(at in c("hsq", "sample_size"))
attr(result, at) <- attr(out_reml, at)
expect_equal(out_reml, result, tol=1e-6)
})
test_that("scan1 works with hsq specified", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
kinship_loco <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
# no covariates, single kinship
out <- scan1(probs, iron$pheno, kinship)
expect_equal(scan1(probs, iron$pheno, kinship, hsq=attr(out, "hsq")), out)
# no covariates, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco)
expect_equal(scan1(probs, iron$pheno, kinship_loco, hsq=attr(out, "hsq")), out)
# X covariates, single kinship
out <- scan1(probs, iron$pheno, kinship, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship, Xcovar=Xc, hsq=attr(out, "hsq")), out)
# X covariates, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship_loco, Xcovar=Xc, hsq=attr(out, "hsq")), out)
# covariates, single kinship
out <- scan1(probs, iron$pheno, kinship, addcovar=X)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, hsq=attr(out, "hsq")), out)
# covariates, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco, addcovar=X)
expect_equal(scan1(probs, iron$pheno, kinship_loco, addcovar=X, hsq=attr(out, "hsq")), out)
# covariates + X covar, single kinship
out <- scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc, hsq=attr(out, "hsq")), out)
# covariates + X covar, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco, addcovar=X, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship_loco, addcovar=X, Xcovar=Xc, hsq=attr(out, "hsq")), out)
})
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("LMM genome scan by scan1 with kinship matrix")
test_that("scan1 with kinship with intercross, vs ported lmmlite code", {
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[1:10,c(1,3)]
map <- insert_pseudomarkers(iron$gmap, step=5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
out_reml <- scan1(probs, iron$pheno, kinship, reml=TRUE)
out_ml <- scan1(probs, iron$pheno, kinship, reml=FALSE)
# "by hand" calculation
y <- iron$pheno
X <- cbind(rep(1, nrow(iron$pheno)))
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% X
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
byhand1_reml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[1,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tol=1e-7)
expect_equal(as.numeric(attr(out_ml, "hsq")[1,]),
c(byhand1_ml$hsq, byhand2_ml$hsq), tol=1e-7)
# compare chromosome 1 LOD scores
d <- dim(probs[[1]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(X, probs[[1]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[1:d,1], lod_reml1, tol=1e-7)
expect_equal(out_reml[1:d,2], lod_reml2, tol=1e-7)
expect_equal(out_ml[1:d,1], lod_ml1, tol=1e-7)
expect_equal(out_ml[1:d,2], lod_ml2, tol=1e-7)
})
test_that("scan1 with intercross with X covariates for null", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
Xc <- get_x_covar(iron)
out_reml <- scan1(probs, iron$pheno, kinship, Xcovar=Xc, reml=TRUE)
out_ml <- scan1(probs, iron$pheno, kinship, Xcovar=Xc, reml=FALSE)
# "by hand" calculation
y <- iron$pheno
X <- cbind(rep(1, nrow(iron$pheno)))
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% X
Xcp <- Ke$vectors %*% Xc
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
byhand1_reml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp), reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp), reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp), reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp), reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[2,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tolerance=1e-6)
expect_equal(as.numeric(attr(out_ml, "hsq")[2,]),
c(byhand1_ml$hsq, byhand2_ml$hsq))
# compare chromosome X LOD scores
d <- dim(probs[["X"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, probs[["X"]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
index <- nrow(out_reml) - rev(1:d) + 1
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1, tolerance=1e-6)
expect_equal(out_reml[index,2], lod_reml2, tolerance=1e-6)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
})
test_that("scan1 with kinship with intercross with an additive covariate", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc, reml=TRUE, tol=1e-12)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc, reml=FALSE, tol=1e-12)
# "by hand" calculation
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% cbind(1, X)
Xcp <- Ke$vectors %*% Xc
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
# autosome null
byhand1A_reml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2A_reml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1A_ml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2A_ml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[1,]),
c(byhand1A_reml$hsq, byhand2A_reml$hsq))
expect_equal(as.numeric(attr(out_ml, "hsq")[1,]),
c(byhand1A_ml$hsq, byhand2A_ml$hsq))
# X chr null
byhand1X_reml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand2X_reml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand1X_ml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
byhand2X_ml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[2,]),
c(byhand1X_reml$hsq, byhand2X_reml$hsq), tolerance=1e-6)
expect_equal(as.numeric(attr(out_ml, "hsq")[2,]),
c(byhand1X_ml$hsq, byhand2X_ml$hsq))
# compare chromosome 2 LOD scores
d <- dim(probs[["2"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, X, probs[["2"]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1A_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2A_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1A_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2A_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1A_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2A_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1A_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2A_ml$loglik)/log(10)
index <- dim(probs[["1"]])[3] + 1:dim(probs[["2"]])[3]
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1)
expect_equal(out_reml[index,2], lod_reml2)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
# compare chromosome X LOD scores
d <- dim(probs[["X"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, probs[["X"]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1X_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2X_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1X_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2X_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1X_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2X_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1X_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2X_ml$loglik)/log(10)
index <- nrow(out_reml) - rev(1:d) + 1
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
## FIX_ME
## REML not yet working on X chromosome, when (X, probs) is not full rank
# expect_equal(out_reml[index,1], lod_reml1)
# expect_equal(out_reml[index,2], lod_reml2)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
})
test_that("scan1 with kinship with intercross with an interactive covariate", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12, intcovar_method="lowmem")
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12, intcovar_method="lowmem")
out_reml_himem <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12, intcovar_method="highmem")
out_ml_himem <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12, intcovar_method="highmem")
# same result using "highmem" and "lowmem" methods
expect_equal(out_reml_himem, out_reml)
expect_equal(out_ml_himem, out_ml)
# "by hand" calculation
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
yp <- Ke$vectors %*% y
Xp <- Ke$vectors %*% cbind(1, X)
Xcp <- Ke$vectors %*% Xc
# double the eigenvalues (== kinship matrix * 2)
Ke$values <- Ke$values*2
# autosome null (same as w/o interactive covariate)
byhand1A_reml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2A_reml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1A_ml <- Rcpp_fitLMM(Ke$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2A_ml <- Rcpp_fitLMM(Ke$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[1,]),
c(byhand1A_reml$hsq, byhand2A_reml$hsq))
expect_equal(as.numeric(attr(out_ml, "hsq")[1,]),
c(byhand1A_ml$hsq, byhand2A_ml$hsq))
# X chr null (same as w/o interactive covariate)
byhand1X_reml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand2X_reml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=TRUE, tol=1e-12)
byhand1X_ml <- Rcpp_fitLMM(Ke$values, yp[,1], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
byhand2X_ml <- Rcpp_fitLMM(Ke$values, yp[,2], cbind(Xp, Xcp[,-1]), reml=FALSE, tol=1e-12)
# hsq the same?
expect_equal(as.numeric(attr(out_reml, "hsq")[2,]),
c(byhand1X_reml$hsq, byhand2X_reml$hsq), tolerance=1e-6)
expect_equal(as.numeric(attr(out_ml, "hsq")[2,]),
c(byhand1X_ml$hsq, byhand2X_ml$hsq))
# compare chromosome 4 LOD scores
npos <- sapply(probs, function(a) dim(a)[[3]])
d <- npos["4"]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, X, probs[["4"]][,-1,i], probs[["4"]][,-1,i]*X)
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1A_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2A_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1A_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2A_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1A_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2A_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1A_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2A_ml$loglik)/log(10)
index <- sum(npos[1:3]) + 1:npos[4]
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1)
expect_equal(out_reml[index,2], lod_reml2)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
# compare chromosome X LOD scores
d <- dim(probs[["X"]])[3]
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, 3)
for(i in 1:d) {
Xp <- Ke$vectors %*% cbind(1, X, probs[["X"]][,-1,i], probs[["X"]][,-1,i]*X)
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1X_reml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2X_reml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1X_ml$hsq, Ke$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2X_ml$hsq, Ke$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1X_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2X_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1X_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2X_ml$loglik)/log(10)
index <- nrow(out_reml) - rev(1:d) + 1
## FIX ME
## Not yet working on X chromosome, when (X, probs) is not full rank
# expect_equal(out_reml[index,1], lod_reml1)
# expect_equal(out_reml[index,2], lod_reml2)
# expect_equal(out_ml[index,1], lod_ml1)
# expect_equal(out_ml[index,2], lod_ml2)
})
test_that("scan1 with kinship works with LOCO, additive covariates", {
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)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
# double the eigenvalues (== kinship matrix * 2)
Ke <- lapply(Ke, function(a) { a$values <- 2*a$values; a})
# compare chromosomes 1, 6, 9, 18
chrs <- paste(c(1,6,9,18))
npos <- sapply(probs, function(a) dim(a)[[3]])
for(chr in chrs) {
nchr <- which(names(npos) == chr)
d <- npos[chr]
yp <- Ke[[chr]]$vectors %*% y
Xp <- Ke[[chr]]$vectors %*% cbind(1, X)
# autosome null
byhand1_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[nchr,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tolerance=1e-5)
expect_equal(as.numeric(attr(out_ml, "hsq")[nchr,]),
c(byhand1_ml$hsq, byhand2_ml$hsq), tolerance=1e-6)
# chromosome scan
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke[[chr]]$vectors %*% cbind(1, X, probs[[chr]][,-1,i])
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
if(nchr > 1) index <- sum(npos[1:(nchr-1)]) + 1:d
else index <- 1:d
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1, tolerance=1e-6)
expect_equal(out_reml[index,2], lod_reml2, tolerance=1e-5)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
}
})
test_that("scan1 with kinship works with LOCO, interactive covariates", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
y <- iron$pheno
Ke <- decomp_kinship(kinship) # eigen decomp
# double the eigenvalues (== kinship matrix * 2)
Ke <- lapply(Ke, function(a) { a$values <- 2*a$values; a})
# compare chromosomes 1, 6, 9, 18
chrs <- paste(c(1,6,9,18))
npos <- sapply(probs, function(a) dim(a)[[3]])
for(chr in chrs) {
nchr <- which(names(npos) == chr)
d <- npos[chr]
yp <- Ke[[chr]]$vectors %*% y
Xp <- Ke[[chr]]$vectors %*% cbind(1, X)
# autosome null (same as w/o interactive covariate)
byhand1_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=TRUE, tol=1e-12)
byhand2_reml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=TRUE, tol=1e-12)
byhand1_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,1], Xp, reml=FALSE, tol=1e-12)
byhand2_ml <- Rcpp_fitLMM(Ke[[chr]]$values, yp[,2], Xp, reml=FALSE, tol=1e-12)
expect_equal(as.numeric(attr(out_reml, "hsq")[nchr,]),
c(byhand1_reml$hsq, byhand2_reml$hsq), tolerance=1e-5)
expect_equal(as.numeric(attr(out_ml, "hsq")[nchr,]),
c(byhand1_ml$hsq, byhand2_ml$hsq), tolerance=1e-6)
# chromosome scan
loglik_reml1 <- loglik_reml2 <-
loglik_ml1 <- loglik_ml2 <- rep(NA, d)
for(i in 1:d) {
Xp <- Ke[[chr]]$vectors %*% cbind(1, X, probs[[chr]][,-1,i], probs[[chr]][,-1,i]*X)
# calculate likelihoods using plain ML (not the residual log likelihood)
loglik_reml1[i] <- Rcpp_calcLL(byhand1_reml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_reml2[i] <- Rcpp_calcLL(byhand2_reml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
loglik_ml1[i] <- Rcpp_calcLL(byhand1_ml$hsq, Ke[[chr]]$values, yp[,1], Xp, reml=FALSE)
loglik_ml2[i] <- Rcpp_calcLL(byhand2_ml$hsq, Ke[[chr]]$values, yp[,2], Xp, reml=FALSE)
}
lod_reml1 <- (loglik_reml1 - byhand1_reml$loglik)/log(10)
lod_reml2 <- (loglik_reml2 - byhand2_reml$loglik)/log(10)
lod_ml1 <- (loglik_ml1 - byhand1_ml$loglik)/log(10)
lod_ml2 <- (loglik_ml2 - byhand2_ml$loglik)/log(10)
if(nchr > 1) index <- sum(npos[1:(nchr-1)]) + 1:d
else index <- 1:d
out_reml <- unclass(out_reml)
out_ml <- unclass(out_ml)
dimnames(out_reml) <- dimnames(out_ml) <- NULL
expect_equal(out_reml[index,1], lod_reml1)
expect_equal(out_reml[index,2], lod_reml2, tolerance=1e-6)
expect_equal(out_ml[index,1], lod_ml1)
expect_equal(out_ml[index,2], lod_ml2)
}
})
test_that("scan1 with kinship works with multicore", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_reml_4core <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12, cores=2)
expect_equal(out_reml, out_reml_4core)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
out_ml_4core <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12, cores=2)
expect_equal(out_ml, out_ml_4core)
})
test_that("scan1 with kinship LOD results invariant to change in scale to pheno and covar", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12)
out_reml_scale <- scan1(probs, iron$pheno/100, kinship, addcovar=X*2, intcovar=X*2,
Xcovar=Xc*2, reml=TRUE, tol=1e-12)
expect_equal(out_reml, out_reml_scale, tol=1e-6)
out_ml <- scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=FALSE, tol=1e-12)
out_ml_scale <- scan1(probs, iron$pheno/100, kinship, addcovar=X*4, intcovar=X*4,
Xcovar=Xc*4, reml=FALSE, tol=1e-12)
expect_equal(out_ml, out_ml_scale, tol=1e-6)
})
test_that("scan1 deals with mismatching 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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
ind <- c(1:50, 101:150)
subK <- lapply(kinship, "[", ind, ind)
expected <- scan1(probs[ind,], iron$pheno[ind,,drop=FALSE], subK, addcovar=X[ind], intcovar=X[ind],
Xcovar=Xc[ind,], reml=TRUE, tol=1e-12)
expect_equal(scan1(probs[ind,], iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno[ind,], kinship, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, subK, addcovar=X, intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X[ind], intcovar=X,
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X[ind],
Xcovar=Xc, reml=TRUE, tol=1e-12), expected)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, intcovar=X,
Xcovar=Xc[ind,], reml=TRUE, tol=1e-12), expected)
})
test_that("scan1 with weights and kinship", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
chr <- c("7", "12", "X")
probs <- probs[,chr]
kinship <- kinship[chr]
RNGkind("Mersenne-Twister") # make sure we're using the standard RNG
set.seed(28915967)
weights <- stats::setNames(sample(1:10, n_ind(iron), replace=TRUE), ind_ids(iron))
out_reml <- scan1(probs, iron$pheno, kinship, addcovar=X,
Xcovar=Xc, weights=weights, reml=TRUE, tol=1e-12)
### results via regress library
# library(regress)
# hsq <- matrix(nrow=3, ncol=2)
# dimnames(hsq) <- list(c("7", "12", "X"), c("liver", "spleen"))
# for(i in 1:2) {
# for(j in 1:3) {
# k <- kinship[[j]]*2
# if(j==3) co <- Xc else co <- X
# out_regress <- regress(iron$pheno[,i] ~ co, ~ k + diag(1/weights), identity=FALSE, tol=1e-8)
# sig <- out_regress$sigma
# hsq[j,i] <- sig[1]/sum(sig)
# }
# }
hsq <- structure(c(0.281878438466624, 0.353558919335662, 0.354315217715963,
0.20945078570429, 0.217671918799751, 0.201682575651958),
.Dim = 3:2, .Dimnames = list(c("7", "12", "X"), c("liver", "spleen")))
expect_equal(attr(out_reml, "hsq"), hsq, tol=1e-5)
### calculate lod scors from scan
# hsq <- attr(out_reml, "hsq")
# result <- vector("list", 6)
# for(i in 1:2) {
# for(j in 1:3) {
# # cholesky decomp of kinhip matrix
# v <- hsq[j,i]*2*kinship[[j]] + (1-hsq[j,i])*diag(1/weights)
# d <- solve(t(chol( v )))
# # pre-multiply phenotype and covar by t(d)
# y <- d %*% iron$pheno[,i]
# if(j==3) co <- Xc else co <- X
# co <- d %*% cbind(1, co)
#
# # log lik under null
# r <- lm(y ~ -1 + co)$resid
# llik0 <- sum(dnorm(r, 0, sqrt(mean(r^2)), log=TRUE))
#
# llik1 <- apply(probs[[j]], 3, function(z) {
# thisX <- cbind(co, d %*% z[,-1,drop=FALSE])
# r <- lm(y ~ -1 + thisX)$resid
# sum(dnorm(r, 0, sqrt(mean(r^2)), log=TRUE)) })
#
# result[[(i-1)*3+j]] <- (llik1 - llik0)/log(10)
#
# }
# }
# result <- rbind(cbind(result[[1]], result[[4]]),
# cbind(result[[2]], result[[5]]),
# cbind(result[[3]], result[[6]]))
# dimnames(result) <- dimnames(out_reml)
# class(result) <- class(out_reml)
result <- structure(c(0.765379546370331, 1.33697830149196, 2.05141217249269,
2.82748389184348, 3.57825928288821, 4.11482576884623,
4.18208695088672, 4.5010472223741, 4.77640847444171,
4.98121670508953, 5.07224467963965, 5.28435499128555,
5.10772158608095, 3.83479814287182, 3.89931966500793,
4.24832760135268, 4.18515409936768, 4.61791967696627,
4.70308738505572, 4.58758777783123, 4.88425659502959,
5.39302146958257, 5.86125538183219, 6.27973204425698,
6.62245498126028, 6.84349072896544, 6.91418168584823,
0.492648140597935, 0.483021280268028, 0.465046322957041,
0.435712712087884, 0.391979816567206, 0.332013881699469,
0.257330299826467, 0.174985513110401, 0.0975448556296765,
0.0390292238934769, 0.0087954970835623, 0.0081071126160765,
0.0315063370227757, 0.0706853835370534, 0.117761588763882,
0.166806834360401, 0.180329873353274, 0.146302909441694,
0.161137797908763, 0.192675853649471, 0.243358929100302,
0.313640413846475, 0.401241925760445, 0.501210880937258,
0.606957139785506, 0.711817493234455, 0.810381904299065,
0.899112752704855, 0.976303366208711, 1.00118015673376,
0.68879116519848, 1.13593060562675, 1.59214070065802,
1.93410070637396, 2.08715532129677, 2.08161469298463,
2.11418090796465, 2.24129850312586, 2.28750006399176,
2.24152015532401, 2.15435105671462, 2.16059659637594,
1.40659065214545, 0.389546568631669, 0.392567989473253,
0.432129729053773, 0.445384879284572, 0.468919719042007,
0.561846602500791, 0.61487493505232, 0.680306211651546,
0.828892839652596, 1.03747954481922, 1.3178288863555,
1.64752359905175, 1.97570834069646, 2.25872692788406,
1.76552836513385, 1.9624748694232, 2.179534029569,
2.40709075366897, 2.62706355895772, 2.81242765604838,
2.93166656075001, 2.95963946310072, 2.89030386259124,
2.74110806468527, 2.54394356128059, 2.33040574729971,
2.12258054254516, 1.93206224497564, 1.76308130810039,
1.61583945469733, 1.57830829020771, 1.91298146944818,
1.89937590028851, 1.86320092206996, 1.79734091344691,
1.69690361501555, 1.56240479342725, 1.4020247631517,
1.23069582035597, 1.06561266027339, 0.920732232061234,
0.803561188099648, 0.715279174632663, 0.690066084958058),
.Dim = c(57L, 2L),
.Dimnames = list(c("D7Mit74", "c7.loc4", "c7.loc6", "c7.loc9",
"c7.loc11", "D7Mit25", "c7.loc14", "c7.loc16",
"c7.loc19", "c7.loc21", "D7Nds5", "c7.loc24",
"c7.loc26", "D7mit30", "c7.loc29", "c7.loc31",
"D7Mit31", "c7.loc34", "c7.loc36", "D7Mit17",
"c7.loc39", "c7.loc41", "c7.loc44", "c7.loc46",
"c7.loc49", "c7.loc51", "D7Mit71", "D12Mit88",
"c12.loc22", "c12.loc25", "c12.loc27", "c12.loc30",
"c12.loc32", "c12.loc35", "c12.loc37", "c12.loc40",
"c12.loc42", "c12.loc45", "c12.loc47", "c12.loc50",
"c12.loc52", "c12.loc55", "c12.loc57", "D12Mit134",
"DXMit16", "cX.loc32", "cX.loc34", "cX.loc37",
"cX.loc40", "cX.loc42", "cX.loc44", "cX.loc47",
"cX.loc50", "cX.loc52", "cX.loc54",
"cX.loc57", "DXMit186"), c("liver", "spleen")),
class = c("scan1", "matrix"))
for(at in c("hsq", "sample_size"))
attr(result, at) <- attr(out_reml, at)
expect_equal(out_reml, result, tol=1e-6)
})
test_that("scan1 works with hsq specified", {
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=2.5)
probs <- calc_genoprob(iron, map, error_prob=0.002)
kinship <- calc_kinship(probs)
kinship_loco <- calc_kinship(probs, "loco")
Xc <- get_x_covar(iron)
X <- match(iron$covar$sex, c("f", "m"))-1
names(X) <- rownames(iron$covar)
# no covariates, single kinship
out <- scan1(probs, iron$pheno, kinship)
expect_equal(scan1(probs, iron$pheno, kinship, hsq=attr(out, "hsq")), out)
# no covariates, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco)
expect_equal(scan1(probs, iron$pheno, kinship_loco, hsq=attr(out, "hsq")), out)
# X covariates, single kinship
out <- scan1(probs, iron$pheno, kinship, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship, Xcovar=Xc, hsq=attr(out, "hsq")), out)
# X covariates, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship_loco, Xcovar=Xc, hsq=attr(out, "hsq")), out)
# covariates, single kinship
out <- scan1(probs, iron$pheno, kinship, addcovar=X)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, hsq=attr(out, "hsq")), out)
# covariates, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco, addcovar=X)
expect_equal(scan1(probs, iron$pheno, kinship_loco, addcovar=X, hsq=attr(out, "hsq")), out)
# covariates + X covar, single kinship
out <- scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship, addcovar=X, Xcovar=Xc, hsq=attr(out, "hsq")), out)
# covariates + X covar, loco kinship
out <- scan1(probs, iron$pheno, kinship_loco, addcovar=X, Xcovar=Xc)
expect_equal(scan1(probs, iron$pheno, kinship_loco, addcovar=X, Xcovar=Xc, hsq=attr(out, "hsq")), out)
})
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.