git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@48927 bc3139a8-67e5-0310-9ffc-ced21a209358
... | ... |
@@ -471,22 +471,16 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){ |
471 | 471 |
W <- W[I, ] |
472 | 472 |
if(any(!is.finite(W))){## | any(!is.finite(V))){ |
473 | 473 |
i <- which(rowSums(!is.finite(W)) > 0) |
474 |
- ##browser() |
|
475 | 474 |
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ") |
476 | 475 |
} |
477 | 476 |
Ns <- tmp.objects[["Ns"]] |
478 | 477 |
Ns <- Ns[I, ] |
479 | 478 |
CHR <- unique(chromosome(object)) |
480 |
- ##batch <- unique(object$batch) |
|
481 | 479 |
batch <- unique(batch(object)) |
482 | 480 |
nuA <- getParam(object, "nuA", batch) |
483 | 481 |
nuB <- getParam(object, "nuB", batch) |
484 |
-## nuA.se <- getParam(object, "nuA.se", batch) |
|
485 |
-## nuB.se <- getParam(object, "nuB.se", batch) |
|
486 | 482 |
phiA <- getParam(object, "phiA", batch) |
487 | 483 |
phiB <- getParam(object, "phiB", batch) |
488 |
-## phiA.se <- getParam(object, "phiA.se", batch) |
|
489 |
-## phiB.se <- getParam(object, "phiB.se", batch) |
|
490 | 484 |
if(CHR==23){ |
491 | 485 |
phiAX <- getParam(object, "phiAX", batch) |
492 | 486 |
phiBX <- getParam(object, "phiBX", batch) |
... | ... |
@@ -494,7 +488,6 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){ |
494 | 488 |
NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05 |
495 | 489 |
if(missing(allele)) stop("must specify allele") |
496 | 490 |
if(CHR == 23){ |
497 |
- ##Design matrix for X chromosome depends on whether there was a sufficient number of AB genotypes |
|
498 | 491 |
if(length(grep("AB", colnames(W))) > 0){ |
499 | 492 |
if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2)) |
500 | 493 |
if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0)) |
... | ... |
@@ -508,27 +501,16 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){ |
508 | 501 |
if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous |
509 | 502 |
betahat <- matrix(NA, 2, nrow(Ystar)) |
510 | 503 |
} |
511 |
- ##How to quickly generate Xstar, Xstar = diag(W) %*% X |
|
512 |
-## Xstar <- apply(W, 1, generateX, X) |
|
513 |
-## IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X)) |
|
514 |
- ##as.numeric(diag(W[1, ]) %*% X) |
|
515 | 504 |
ww <- rep(1, ncol(Ystar)) |
516 | 505 |
II <- which(rowSums(is.nan(Ystar)) == 0) |
517 | 506 |
for(i in II){ |
518 | 507 |
betahat[, i] <- dqrlsWrapper(W[i, ] * X, Ystar[i, ], ww) |
519 |
- ##betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ])) |
|
520 |
- ##ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2) |
|
521 |
- ##ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr)) |
|
522 | 508 |
} |
523 | 509 |
if(allele == "A"){ |
524 | 510 |
nuA[complete] <- betahat[1, ] |
525 | 511 |
phiA[complete] <- betahat[2, ] |
526 |
-## nuA.se[complete] <- ses[1, ] |
|
527 |
-## phiA.se[complete] <- ses[2, ] |
|
528 | 512 |
object <- pr(object, "nuA", batch, nuA) |
529 |
-## object <- pr(object, "nuA.se", batch, nuA.se) |
|
530 | 513 |
object <- pr(object, "phiA", batch, phiA) |
531 |
-## object <- pr(object, "phiA.se", batch, phiA.se) |
|
532 | 514 |
if(CHR == 23){ |
533 | 515 |
phiAX[complete] <- betahat[3, ] |
534 | 516 |
object <- pr(object, "phiAX", batch, phiAX) |
... | ... |
@@ -537,23 +519,13 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){ |
537 | 519 |
if(allele=="B"){ |
538 | 520 |
nuB[complete] <- betahat[1, ] |
539 | 521 |
phiB[complete] <- betahat[2, ] |
540 |
-## nuB.se[complete] <- ses[1, ] |
|
541 |
-## phiB.se[complete] <- ses[2, ] |
|
542 | 522 |
object <- pr(object, "nuB", batch, nuB) |
543 |
-## object <- pr(object, "nuB.se", batch, nuB.se) |
|
544 | 523 |
object <- pr(object, "phiB", batch, phiB) |
545 |
-## object <- pr(object, "phiB.se", batch, phiB.se) |
|
546 | 524 |
if(CHR == 23){ |
547 | 525 |
phiBX[complete] <- betahat[3, ] |
548 | 526 |
object <- pr(object, "phiBX", batch, phiBX) |
549 | 527 |
} |
550 | 528 |
} |
551 |
-## THR.NU.PHI <- cnOptions$THR.NU.PHI |
|
552 |
-## if(THR.NU.PHI){ |
|
553 |
-## verbose <- cnOptions$verbose |
|
554 |
-## if(verbose) message("Thresholding nu and phi") |
|
555 |
-## object <- thresholdModelParams(object, cnOptions) |
|
556 |
-## } |
|
557 | 529 |
return(object) |
558 | 530 |
} |
559 | 531 |
|