R/kinship.on.the.fly.R

Defines functions kinship.on.the.fly

Documented in kinship.on.the.fly

kinship.on.the.fly <-
function(kin.mats, geno, chr1 = NULL, chr2 = NULL, phenoV, covarV = NULL){

	get.g=function(pair = NULL, phenotype, covarV){
		
		if(is.null(pair)){
			pair.name <- "full.geno"
			}else{
			pair.name <- paste(pair, collapse = ",")
			}
		kin.mat.locale <- which(names(kin.mats) == pair.name)
		K <- kin.mats[[kin.mat.locale]]
		
		#if we are correcting the covariate only don't put it in the model
		if(is.null(covarV) || is.null(pair)){
			model = regress(as.vector(phenotype)~1,~K, pos = c(TRUE, TRUE))	
			}else{
			model = regress(as.vector(phenotype)~covarV, ~K, pos = c(TRUE, TRUE))		
			}
				
		#This err.cov is the same as err.cov in Dan's code using estVC
		err.cov = summary(model)$sigma[1]*K+summary(model)$sigma[2]*diag(nrow(K))

		
		eW = eigen(err.cov, symmetric = TRUE)
	    if(min(eW$values) < 0 && abs(min(eW$values)) > sqrt(.Machine$double.eps)){
		      }else{
	      	eW$values[eW$values <= 0] = Inf
	      	} 
	      	err.cov = eW$vector %*% diag(eW$values^-0.5) %*% t(eW$vector)

		new.pheno <- err.cov %*% phenotype
		new.geno <- err.cov %*% geno; #hist(new.geno)
		rownames(new.geno) <- rownames(geno)
		if(!is.null(covarV)){
			new.covar <- err.cov %*% covarV
			}else{
			new.covar <- NULL	
			}
		
		results = list(err.cov, new.pheno, new.geno, new.covar)
		names(results) <- c("err.cov", "corrected.pheno", "corrected.geno", "corrected.covar")
		return(results)
		}

	result <- get.g(pair = c(chr1, chr2), phenotype = phenoV, covarV = covarV)
	
	return(result)
	
}

Try the cape package in your browser

Any scripts or data that you put into this service are public.

cape documentation built on May 2, 2019, 3:27 a.m.