Mercurial > repos > dereeper > mlmm
diff mlmm/source_library/emma.mlmm.r @ 0:6b7107812931 draft
Uploaded
| author | dereeper |
|---|---|
| date | Thu, 02 Jul 2015 05:42:38 -0400 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mlmm/source_library/emma.mlmm.r Thu Jul 02 05:42:38 2015 -0400 @@ -0,0 +1,1274 @@ +emma.kinship <- function(snps, method="additive", use="all") { + n0 <- sum(snps==0,na.rm=TRUE) + nh <- sum(snps==0.5,na.rm=TRUE) + n1 <- sum(snps==1,na.rm=TRUE) + nNA <- sum(is.na(snps)) + + stopifnot(n0+nh+n1+nNA == length(snps)) + + if ( method == "dominant" ) { + flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps)) + snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)] + } + else if ( method == "recessive" ) { + flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps)) + snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)] + } + else if ( ( method == "additive" ) && ( nh > 0 ) ) { + dsnps <- snps + rsnps <- snps + flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps)) + dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)] + flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps)) + rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)] + snps <- rbind(dsnps,rsnps) + } + + if ( use == "all" ) { + mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps)) + snps[is.na(snps)] <- mafs[is.na(snps)] + } + else if ( use == "complete.obs" ) { + snps <- snps[rowSums(is.na(snps))==0,] + } + + n <- ncol(snps) + K <- matrix(nrow=n,ncol=n) + diag(K) <- 1 + + for(i in 2:n) { + for(j in 1:(i-1)) { + x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j]) + K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x)) + K[j,i] <- K[i,j] + } + } + return(K) +} + +emma.eigen.L <- function(Z,K,complete=TRUE) { + if ( is.null(Z) ) { + return(emma.eigen.L.wo.Z(K)) + } + else { + return(emma.eigen.L.w.Z(Z,K,complete)) + } +} + +emma.eigen.L.wo.Z <- function(K) { + eig <- eigen(K,symmetric=TRUE) + return(list(values=eig$values,vectors=eig$vectors)) +} + +emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) { + if ( complete == FALSE ) { + vids <- colSums(Z)>0 + Z <- Z[,vids] + K <- K[vids,vids] + } + eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE) + return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE))) +} + +emma.eigen.R <- function(Z,K,X,complete=TRUE) { + if ( ncol(X) == 0 ) { + return(emma.eigen.L(Z,K)) + } + else if ( is.null(Z) ) { + return(emma.eigen.R.wo.Z(K,X)) + } + else { + return(emma.eigen.R.w.Z(Z,K,X,complete)) + } +} + +emma.eigen.R.wo.Z <- function(K, X) { + n <- nrow(X) + q <- ncol(X) + S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X) + eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE) + stopifnot(!is.complex(eig$values)) + return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)])) +} + +emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) { + if ( complete == FALSE ) { + vids <- colSums(Z) > 0 + Z <- Z[,vids] + K <- K[vids,vids] + } + n <- nrow(Z) + t <- ncol(Z) + q <- ncol(X) + + SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z) + eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE) + if ( is.complex(eig$values) ) { + eig$values <- Re(eig$values) + eig$vectors <- Re(eig$vectors) + } + qr.X <- qr.Q(qr(X)) + return(list(values=eig$values[1:(t-q)], + vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)), + complete=TRUE)[,c(1:(t-q),(t+1):n)])) +} + +emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) { + n <- length(xi) + delta <- exp(logdelta) + return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) ) +} + +emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) { + t <- length(xi.1) + delta <- exp(logdelta) +# stopifnot(length(lambda) == length(etas.1)) + return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) ) +} + +emma.delta.ML.dLL.wo.Z <- function(logdelta, lambda, etas, xi) { + n <- length(xi) + delta <- exp(logdelta) + etasq <- etas*etas + ldelta <- lambda+delta + return( 0.5*(n*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/(xi+delta))) ) +} + +emma.delta.ML.dLL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) { + t <- length(xi.1) + delta <- exp(logdelta) + etasq <- etas.1*etas.1 + ldelta <- lambda+delta + return( 0.5*(n*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/(xi.1+delta))+(n-t)/delta) ) ) +} + +emma.delta.REML.LL.wo.Z <- function(logdelta, lambda, etas) { + nq <- length(etas) + delta <- exp(logdelta) + return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas*etas/(lambda+delta))))-sum(log(lambda+delta))) ) +} + +emma.delta.REML.LL.w.Z <- function(logdelta, lambda, etas.1, n, t, etas.2.sq ) { + tq <- length(etas.1) + nq <- n - t + tq + delta <- exp(logdelta) + return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(lambda+delta))+(n-t)*logdelta)) ) +} + +emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) { + nq <- length(etas) + delta <- exp(logdelta) + etasq <- etas*etas + ldelta <- lambda+delta + return( 0.5*(nq*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/ldelta)) ) +} + +emma.delta.REML.dLL.w.Z <- function(logdelta, lambda, etas.1, n, t1, etas.2.sq ) { + t <- t1 + tq <- length(etas.1) + nq <- n - t + tq + delta <- exp(logdelta) + etasq <- etas.1*etas.1 + ldelta <- lambda+delta + return( 0.5*(nq*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/ldelta)+(n-t)/delta)) ) +} + +emma.MLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10, + esp=1e-10, eig.L = NULL, eig.R = NULL) +{ + n <- length(y) + t <- nrow(K) + q <- ncol(X) + +# stopifnot(nrow(K) == t) + stopifnot(ncol(K) == t) + stopifnot(nrow(X) == n) + + if ( det(crossprod(X,X)) == 0 ) { + warning("X is singular") + return (list(ML=0,delta=0,ve=0,vg=0)) + } + + if ( is.null(Z) ) { + if ( is.null(eig.L) ) { + eig.L <- emma.eigen.L.wo.Z(K) + } + if ( is.null(eig.R) ) { + eig.R <- emma.eigen.R.wo.Z(K,X) + } + etas <- crossprod(eig.R$vectors,y) + + + logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim + m <- length(logdelta) + delta <- exp(logdelta) + Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE) + Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE) + Etasq <- matrix(etas*etas,n-q,m) + LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Xis))) + dLL <- 0.5*delta*(n*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Xis)) + + optlogdelta <- vector(length=0) + optLL <- vector(length=0) + if ( dLL[1] < esp ) { + optlogdelta <- append(optlogdelta, llim) + optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.R$values,etas,eig.L$values)) + } + if ( dLL[m-1] > 0-esp ) { + optlogdelta <- append(optlogdelta, ulim) + optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.R$values,etas,eig.L$values)) + } + + for( i in 1:(m-1) ) + { + if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + { + r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas, xi=eig.L$values) + optlogdelta <- append(optlogdelta, r$root) + optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.R$values, etas, eig.L$values)) + } + } +# optdelta <- exp(optlogdelta) + } + else { + if ( is.null(eig.L) ) { + eig.L <- emma.eigen.L.w.Z(Z,K) + } + if ( is.null(eig.R) ) { + eig.R <- emma.eigen.R.w.Z(Z,K,X) + } + etas <- crossprod(eig.R$vectors,y) + etas.1 <- etas[1:(t-q)] + etas.2 <- etas[(t-q+1):(n-q)] + etas.2.sq <- sum(etas.2*etas.2) + + logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim + + m <- length(logdelta) + delta <- exp(logdelta) + Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE) + Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE) + Etasq <- matrix(etas.1*etas.1,t-q,m) + #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas)) + dLL <- 0.5*delta*(n*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta)) + + optlogdelta <- vector(length=0) + optLL <- vector(length=0) + if ( dLL[1] < esp ) { + optlogdelta <- append(optlogdelta, llim) + optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq)) + } + if ( dLL[m-1] > 0-esp ) { + optlogdelta <- append(optlogdelta, ulim) + optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq)) + } + + for( i in 1:(m-1) ) + { + if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + { + r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq ) + optlogdelta <- append(optlogdelta, r$root) + optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.R$values, etas.1, eig.L$values, n, etas.2.sq )) + } + } +# optdelta <- exp(optlogdelta) + } + + maxdelta <- exp(optlogdelta[which.max(optLL)]) + maxLL <- max(optLL) + if ( is.null(Z) ) { + maxva <- sum(etas*etas/(eig.R$values+maxdelta))/n + } + else { + maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/n + } + maxve <- maxva*maxdelta + + return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva)) +} + +emma.MLE.noX <- function(y, K, Z=NULL, ngrids=100, llim=-10, ulim=10, + esp=1e-10, eig.L = NULL) +{ + n <- length(y) + t <- nrow(K) + +# stopifnot(nrow(K) == t) + stopifnot(ncol(K) == t) + + if ( is.null(Z) ) { + if ( is.null(eig.L) ) { + eig.L <- emma.eigen.L.wo.Z(K) + } + etas <- crossprod(eig.L$vectors,y) + + logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim + m <- length(logdelta) + delta <- exp(logdelta) + Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE) + Etasq <- matrix(etas*etas,n,m) + LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Xis)))-colSums(log(Xis))) + dLL <- 0.5*delta*(n*colSums(Etasq/(Xis*Xis))/colSums(Etasq/Xis)-colSums(1/Xis)) + + optlogdelta <- vector(length=0) + optLL <- vector(length=0) + #print(dLL) + if ( dLL[1] < esp ) { + optlogdelta <- append(optlogdelta, llim) + optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.L$values,etas,eig.L$values)) + } + if ( dLL[m-1] > 0-esp ) { + optlogdelta <- append(optlogdelta, ulim) + optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.L$values,etas,eig.L$values)) + } + + for( i in 1:(m-1) ) + { + #if ( ( dLL[i]*dLL[i+1] < 0 ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + { + r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas=etas, xi=eig.L$values) + optlogdelta <- append(optlogdelta, r$root) + optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.L$values, etas, eig.L$values)) + } + } +# optdelta <- exp(optlogdelta) + } + else { + if ( is.null(eig.L) ) { + eig.L <- emma.eigen.L.w.Z(Z,K) + } + etas <- crossprod(eig.L$vectors,y) + etas.1 <- etas[1:t] + etas.2 <- etas[(t+1):n] + etas.2.sq <- sum(etas.2*etas.2) + + logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim + + m <- length(logdelta) + delta <- exp(logdelta) + Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE) + Etasq <- matrix(etas.1*etas.1,t,m) + #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas)) + dLL <- 0.5*delta*(n*(colSums(Etasq/(Xis*Xis))+etas.2.sq/(delta*delta))/(colSums(Etasq/Xis)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta)) + + optlogdelta <- vector(length=0) + optLL <- vector(length=0) + if ( dLL[1] < esp ) { + optlogdelta <- append(optlogdelta, llim) + optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq)) + } + if ( dLL[m-1] > 0-esp ) { + optlogdelta <- append(optlogdelta, ulim) + optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq)) + } + + for( i in 1:(m-1) ) + { + if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + { + r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq ) + optlogdelta <- append(optlogdelta, r$root) + optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.L$values, etas.1, eig.L$values, n, etas.2.sq )) + } + } +# optdelta <- exp(optlogdelta) + } + + maxdelta <- exp(optlogdelta[which.max(optLL)]) + maxLL <- max(optLL) + if ( is.null(Z) ) { + maxva <- sum(etas*etas/(eig.L$values+maxdelta))/n + } + else { + maxva <- (sum(etas.1*etas.1/(eig.L$values+maxdelta))+etas.2.sq/maxdelta)/n + } + maxve <- maxva*maxdelta + + return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva)) +} + +emma.REMLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10, + esp=1e-10, eig.L = NULL, eig.R = NULL) { + n <- length(y) + t <- nrow(K) + q <- ncol(X) + +# stopifnot(nrow(K) == t) + stopifnot(ncol(K) == t) + stopifnot(nrow(X) == n) + + if ( det(crossprod(X,X)) == 0 ) { + warning("X is singular") + return (list(REML=0,delta=0,ve=0,vg=0)) + } + + if ( is.null(Z) ) { + if ( is.null(eig.R) ) { + eig.R <- emma.eigen.R.wo.Z(K,X) + } + etas <- crossprod(eig.R$vectors,y) + + logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim + m <- length(logdelta) + delta <- exp(logdelta) + Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE) + Etasq <- matrix(etas*etas,n-q,m) + LL <- 0.5*((n-q)*(log((n-q)/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Lambdas))) + dLL <- 0.5*delta*((n-q)*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Lambdas)) + + optlogdelta <- vector(length=0) + optLL <- vector(length=0) + if ( dLL[1] < esp ) { + optlogdelta <- append(optlogdelta, llim) + optLL <- append(optLL, emma.delta.REML.LL.wo.Z(llim,eig.R$values,etas)) + } + if ( dLL[m-1] > 0-esp ) { + optlogdelta <- append(optlogdelta, ulim) + optLL <- append(optLL, emma.delta.REML.LL.wo.Z(ulim,eig.R$values,etas)) + } + + for( i in 1:(m-1) ) + { + if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + { + r <- uniroot(emma.delta.REML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas) + optlogdelta <- append(optlogdelta, r$root) + optLL <- append(optLL, emma.delta.REML.LL.wo.Z(r$root,eig.R$values, etas)) + } + } +# optdelta <- exp(optlogdelta) + } + else { + if ( is.null(eig.R) ) { + eig.R <- emma.eigen.R.w.Z(Z,K,X) + } + etas <- crossprod(eig.R$vectors,y) + etas.1 <- etas[1:(t-q)] + etas.2 <- etas[(t-q+1):(n-q)] + etas.2.sq <- sum(etas.2*etas.2) + + logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim + m <- length(logdelta) + delta <- exp(logdelta) + Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE) + Etasq <- matrix(etas.1*etas.1,t-q,m) + dLL <- 0.5*delta*((n-q)*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Lambdas)+(n-t)/delta)) + + optlogdelta <- vector(length=0) + optLL <- vector(length=0) + if ( dLL[1] < esp ) { + optlogdelta <- append(optlogdelta, llim) + optLL <- append(optLL, emma.delta.REML.LL.w.Z(llim,eig.R$values,etas.1,n,t,etas.2.sq)) + } + if ( dLL[m-1] > 0-esp ) { + optlogdelta <- append(optlogdelta, ulim) + optLL <- append(optLL, emma.delta.REML.LL.w.Z(ulim,eig.R$values,etas.1,n,t,etas.2.sq)) + } + + for( i in 1:(m-1) ) + { + if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) + { + r <- uniroot(emma.delta.REML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, n=n, t1=t, etas.2.sq = etas.2.sq ) + optlogdelta <- append(optlogdelta, r$root) + optLL <- append(optLL, emma.delta.REML.LL.w.Z(r$root,eig.R$values, etas.1, n, t, etas.2.sq )) + } + } +# optdelta <- exp(optlogdelta) + } + + maxdelta <- exp(optlogdelta[which.max(optLL)]) + maxLL <- max(optLL) + if ( is.null(Z) ) { + maxva <- sum(etas*etas/(eig.R$values+maxdelta))/(n-q) + } + else { + maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/(n-q) + } + maxve <- maxva*maxdelta + + return (list(REML=maxLL,delta=maxdelta,ve=maxve,vg=maxva)) +} + +emma.ML.LRT <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) { + if ( is.null(dim(ys)) || ncol(ys) == 1 ) { + ys <- matrix(ys,1,length(ys)) + } + if ( is.null(dim(xs)) || ncol(xs) == 1 ) { + xs <- matrix(xs,1,length(xs)) + } + if ( is.null(X0) ) { + X0 <- matrix(1,ncol(ys),1) + } + + g <- nrow(ys) + n <- ncol(ys) + m <- nrow(xs) + t <- ncol(xs) + q0 <- ncol(X0) + q1 <- q0 + 1 + + if ( !ponly ) { + ML1s <- matrix(nrow=m,ncol=g) + ML0s <- matrix(nrow=m,ncol=g) + vgs <- matrix(nrow=m,ncol=g) + ves <- matrix(nrow=m,ncol=g) + } + stats <- matrix(nrow=m,ncol=g) + ps <- matrix(nrow=m,ncol=g) + ML0 <- vector(length=g) + + stopifnot(nrow(K) == t) + stopifnot(ncol(K) == t) + stopifnot(nrow(X0) == n) + + if ( sum(is.na(ys)) == 0 ) { + eig.L <- emma.eigen.L(Z,K) + eig.R0 <- emma.eigen.R(Z,K,X0) + + for(i in 1:g) { + ML0[i] <- emma.MLE(ys[i,],X0,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R0)$ML + } + + x.prev <- vector(length=0) + + for(i in 1:m) { + vids <- !is.na(xs[i,]) + nv <- sum(vids) + xv <- xs[i,vids] + + if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { + if (!ponly) { + stats[i,] <- rep(NA,g) + vgs[i,] <- rep(NA,g) + ves[i,] <- rep(NA,g) + ML1s[i,] <- rep(NA,g) + ML0s[i,] <- rep(NA,g) + } + ps[i,] = rep(1,g) + } + else if ( identical(x.prev, xv) ) { + if ( !ponly ) { + stats[i,] <- stats[i-1,] + vgs[i,] <- vgs[i-1,] + ves[i,] <- ves[i-1,] + ML1s[i,] <- ML1s[i-1,] + ML0s[i,] <- ML0s[i-1,] + } + ps[i,] <- ps[i-1,] + } + else { + if ( is.null(Z) ) { + X <- cbind(X0[vids,,drop=FALSE],xs[i,vids]) + eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X) + } + else { + vrows <- as.logical(rowSums(Z[,vids])) + nr <- sum(vrows) + X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids]%*%t(xs[i,vids,drop=FALSE])) + eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X) + } + + for(j in 1:g) { + if ( nv == t ) { + MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1) +# MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1) + if (!ponly) { + ML1s[i,j] <- MLE$ML + vgs[i,j] <- MLE$vg + ves[i,j] <- MLE$ve + } + stats[i,j] <- 2*(MLE$ML-ML0[j]) + + } + else { + if ( is.null(Z) ) { + eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids]) + MLE0 <- emma.MLE(ys[j,vids],X0[vids,,drop=FALSE],K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0) + MLE1 <- emma.MLE(ys[j,vids],X,K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0) + } + else { + if ( nr == n ) { + MLE1 <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L) + } + else { + eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids],K[vids,vids]) + MLE0 <- emma.MLE(ys[j,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0) + MLE1 <- emma.MLE(ys[j,vrows],X,K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0) + } + } + if (!ponly) { + ML1s[i,j] <- MLE1$ML + ML0s[i,j] <- MLE0$ML + vgs[i,j] <- MLE1$vg + ves[i,j] <- MLE1$ve + } + stats[i,j] <- 2*(MLE1$ML-MLE0$ML) + } + } + if ( ( nv == t ) && ( !ponly ) ) { + ML0s[i,] <- ML0 + } + ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE) + } + } + } + else { + eig.L <- emma.eigen.L(Z,K) + eig.R0 <- emma.eigen.R(Z,K,X0) + + for(i in 1:g) { + vrows <- !is.na(ys[i,]) + if ( is.null(Z) ) { + ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vrows,vrows],NULL,ngrids,llim,ulim,esp)$ML + } + else { + vids <- colSums(Z[vrows,]>0) + + ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp)$ML + } + } + + x.prev <- vector(length=0) + + for(i in 1:m) { + vids <- !is.na(xs[i,]) + nv <- sum(vids) + xv <- xs[i,vids] + + if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { + if (!ponly) { + stats[i,] <- rep(NA,g) + vgs[i,] <- rep(NA,g) + ves[i,] <- rep(NA,g) + ML1s[i,] <- rep(NA,g) + ML0s[,i] <- rep(NA,g) + } + ps[i,] = rep(1,g) + } + else if ( identical(x.prev, xv) ) { + if ( !ponly ) { + stats[i,] <- stats[i-1,] + vgs[i,] <- vgs[i-1,] + ves[i,] <- ves[i-1,] + ML1s[i,] <- ML1s[i-1,] + } + ps[i,] = ps[i-1,] + } + else { + if ( is.null(Z) ) { + X <- cbind(X0,xs[i,]) + if ( nv == t ) { + eig.R1 = emma.eigen.R.wo.Z(K,X) + } + } + else { + vrows <- as.logical(rowSums(Z[,vids])) + X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) + if ( nv == t ) { + eig.R1 = emma.eigen.R.w.Z(Z,K,X) + } + } + + for(j in 1:g) { +# print(j) + vrows <- !is.na(ys[j,]) + if ( nv == t ) { + nr <- sum(vrows) + if ( is.null(Z) ) { + if ( nr == n ) { + MLE <- emma.MLE(ys[j,],X,K,NULL,ngrids,llim,ulim,esp,eig.L,eig.R1) + } + else { + MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vrows,vrows],NULL,ngrids,llim,ulim,esp) + } + } + else { + if ( nr == n ) { + MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1) + } + else { + vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) + MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vtids,vtids],Z[vrows,vtids],ngrids,llim,ulim,esp) + } + } + + if (!ponly) { + ML1s[i,j] <- MLE$ML + vgs[i,j] <- MLE$vg + ves[i,j] <- MLE$ve + } + stats[i,j] <- 2*(MLE$ML-ML0[j]) + } + else { + if ( is.null(Z) ) { + vtids <- vrows & vids + eig.L0 <- emma.eigen.L(NULL,K[vtids,vtids]) + MLE0 <- emma.MLE(ys[j,vtids],X0[vtids,,drop=FALSE],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0) + MLE1 <- emma.MLE(ys[j,vtids],X[vtids,],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0) + } + else { + vtids <- as.logical(colSums(Z[vrows,])) & vids + vtrows <- vrows & as.logical(rowSums(Z[,vids])) + eig.L0 <- emma.eigen.L(Z[vtrows,vtids],K[vtids,vtids]) + MLE0 <- emma.MLE(ys[j,vtrows],X0[vtrows,,drop=FALSE],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0) + MLE1 <- emma.MLE(ys[j,vtrows],X[vtrows,],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0) + } + if (!ponly) { + ML1s[i,j] <- MLE1$ML + vgs[i,j] <- MLE1$vg + ves[i,j] <- MLE1$ve + ML0s[i,j] <- MLE0$ML + } + stats[i,j] <- 2*(MLE1$ML-MLE0$ML) + } + } + if ( ( nv == t ) && ( !ponly ) ) { + ML0s[i,] <- ML0 + } + ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE) + } + } + } + if ( ponly ) { + return (ps) + } + else { + return (list(ps=ps,ML1s=ML1s,ML0s=ML0s,stats=stats,vgs=vgs,ves=ves)) + } +} + +emma.test <- function(ys, xs, K, Z=NULL, x0s = NULL, X0 = NULL, dfxs = 1, dfx0s = 1, use.MLE = FALSE, use.LRT = FALSE, ngrids = 100, llim = -10, ulim = 10, esp=1e-10, ponly = FALSE) +{ + stopifnot (dfxs > 0) + + if ( is.null(dim(ys)) || ncol(ys) == 1 ) { + ys <- matrix(ys,1,length(ys)) + } + + if ( is.null(dim(xs)) || ncol(xs) == 1 ) { + xs <- matrix(xs,1,length(xs)) + } + nx <- nrow(xs)/dfxs + + if ( is.null(x0s) ) { + dfx0s = 0 + x0s <- matrix(NA,0,ncol(xs)) + } + # X0 automatically contains intercept. If no intercept is to be used, + # X0 should be matrix(nrow=ncol(ys),ncol=0) + if ( is.null(X0) ) { + X0 <- matrix(1,ncol(ys),1) + } + + stopifnot(Z == NULL) # The case where Z is not null is not implemented + + ny <- nrow(ys) + iy <- ncol(ys) + ix <- ncol(xs) + + stopifnot(nrow(K) == ix) + stopifnot(ncol(K) == ix) + stopifnot(nrow(X0) == iy) + + if ( !ponly ) { + LLs <- matrix(nrow=m,ncol=g) + vgs <- matrix(nrow=m,ncol=g) + ves <- matrix(nrow=m,ncol=g) + } + dfs <- matrix(nrow=m,ncol=g) + stats <- matrix(nrow=m,ncol=g) + ps <- matrix(nrow=m,ncol=g) + + # The case with no missing phenotypes + if ( sum(is.na(ys)) == 0 ) { + if ( ( use.MLE ) || ( !use.LRT ) ) { + eig.L0 <- emma.eigen.L(Z,K) + } + if ( dfx0s == 0 ) { + eig.R0 <- emma.eigen.R(Z,K,X0) + } + x.prev <- NULL + + for(i in 1:ix) { + x1 <- t(xs[(dfxs*(i-1)+1):(dfxs*i),,drop=FALSE]) + if ( dfxs0 == 0 ) { + x0 <- X0 + } + else { + x0 <- cbind(t(x0s[(dfx0s*(i-1)+1):(dfx0s*i),,drop=FALSE]),X0) + } + x <- cbind(x1,x0) + xvids <- rowSums(is.na(x) == 0) + nxv <- sum(xvids) + xv <- x[xvids,,drop=FALSE] + Kv <- K[xvids,xvids,drop=FALSE] + yv <- ys[j,xvids] + + if ( identical(x.prev, xv) ) { + if ( !ponly ) { + vgs[i,] <- vgs[i-1,] + ves[i,] <- ves[i-1,] + dfs[i,] <- dfs[i-1,] + REMLs[i,] <- REMLs[i-1,] + stats[i,] <- stats[i-1,] + } + ps[i,] <- ps[i-1,] + } + else { + eig.R1 = emma.eigen.R.wo.Z(Kv,xv) + + for(j in 1:iy) { + if ( ( use.MLE ) || ( !use.LRT ) ) { + if ( nxv < t ) { + # NOTE: this complexity can be improved by avoiding eigen computation for identical missing patterns + eig.L0v <- emma.eigen.L.wo.Z(Kv) + } + else { + eig.L0v <- eig.L0 + } + } + + if ( use.MLE ) { + MLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1) + stop("Not implemented yet") + } + else { + REMLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1) + if ( use.LRT ) { + stop("Not implemented yet") + } + else { + U <- eig.L0v$vectors * matrix(sqrt(1/(eig.L0v$values+REMLE$delta)),t,t,byrow=TRUE) + dfs[i,j] <- length(eig.R1$values) + yt <- crossprod(U,yv) + xt <- crossprod(U,xv) + ixx <- solve(crossprod(xt,xt)) + beta <- ixx%*%crossprod(xt,yt) + if ( dfxs == 1 ) { + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + } + else { + model.m <- c(rep(1,dfxs),rep(0,ncol(xv)-dfxs)) + stats[i,j] <- + crossprod(crossprod(solve(crossprod(crossprod(iXX,model.m), + model.m)), + model.m*beta),model.m*beta) + + } + if ( !ponly ) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + } + } + } + if ( dfxs == 1 ) { + ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) + } + else { + ps[i,] <- pf(abs(stats[i,]),dfs[i,],lower.tail=FALSE) + } + } + } + } + # The case with missing genotypes - not implemented yet + else { + stop("Not implemented yet") + eig.L <- emma.eigen.L(Z,K) + eig.R0 <- emma.eigen.R(Z,K,X0) + + x.prev <- vector(length=0) + + for(i in 1:m) { + vids <- !is.na(xs[i,]) + nv <- sum(vids) + xv <- xs[i,vids] + + if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { + if (!ponly) { + vgs[i,] <- rep(NA,g) + ves[i,] <- rep(NA,g) + REMLs[i,] <- rep(NA,g) + dfs[i,] <- rep(NA,g) + } + ps[i,] = rep(1,g) + } + else if ( identical(x.prev, xv) ) { + if ( !ponly ) { + stats[i,] <- stats[i-1,] + vgs[i,] <- vgs[i-1,] + ves[i,] <- ves[i-1,] + REMLs[i,] <- REMLs[i-1,] + dfs[i,] <- dfs[i-1,] + } + ps[i,] = ps[i-1,] + } + else { + if ( is.null(Z) ) { + X <- cbind(X0,xs[i,]) + if ( nv == t ) { + eig.R1 = emma.eigen.R.wo.Z(K,X) + } + } + else { + vrows <- as.logical(rowSums(Z[,vids,drop=FALSE])) + X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) + if ( nv == t ) { + eig.R1 = emma.eigen.R.w.Z(Z,K,X) + } + } + + for(j in 1:g) { + vrows <- !is.na(ys[j,]) + if ( nv == t ) { + yv <- ys[j,vrows] + nr <- sum(vrows) + if ( is.null(Z) ) { + if ( nr == n ) { + REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1) + U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) + } + else { + eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE]) + REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) + } + dfs[i,j] <- nr-q1 + } + else { + if ( nr == n ) { + REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1) + U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) + } + else { + vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) + eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) + REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) + } + dfs[i,j] <- nr-q1 + } + + yt <- crossprod(U,yv) + Xt <- crossprod(U,X[vrows,,drop=FALSE]) + iXX <- solve(crossprod(Xt,Xt)) + beta <- iXX%*%crossprod(Xt,yt) + if ( !ponly ) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + } + else { + if ( is.null(Z) ) { + vtids <- vrows & vids + eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE]) + yv <- ys[j,vtids] + nr <- sum(vtids) + REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) + Xt <- crossprod(U,X[vtids,,drop=FALSE]) + dfs[i,j] <- nr-q1 + } + else { + vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids + vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE])) + eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) + yv <- ys[j,vtrows] + nr <- sum(vtrows) + REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) + Xt <- crossprod(U,X[vtrows,,drop=FALSE]) + dfs[i,j] <- nr-q1 + } + yt <- crossprod(U,yv) + iXX <- solve(crossprod(Xt,Xt)) + beta <- iXX%*%crossprod(Xt,yt) + if ( !ponly ) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + + } + } + ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) + } + } + } + if ( ponly ) { + return (ps) + } + else { + return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves)) + } +} + +emma.REML.t <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) { + if ( is.null(dim(ys)) || ncol(ys) == 1 ) { + ys <- matrix(ys,1,length(ys)) + } + if ( is.null(dim(xs)) || ncol(xs) == 1 ) { + xs <- matrix(xs,1,length(xs)) + } + if ( is.null(X0) ) { + X0 <- matrix(1,ncol(ys),1) + } + + g <- nrow(ys) + n <- ncol(ys) + m <- nrow(xs) + t <- ncol(xs) + q0 <- ncol(X0) + q1 <- q0 + 1 + + stopifnot(nrow(K) == t) + stopifnot(ncol(K) == t) + stopifnot(nrow(X0) == n) + + if ( !ponly ) { + REMLs <- matrix(nrow=m,ncol=g) + vgs <- matrix(nrow=m,ncol=g) + ves <- matrix(nrow=m,ncol=g) + } + dfs <- matrix(nrow=m,ncol=g) + stats <- matrix(nrow=m,ncol=g) + ps <- matrix(nrow=m,ncol=g) + + if ( sum(is.na(ys)) == 0 ) { + eig.L <- emma.eigen.L(Z,K) + + x.prev <- vector(length=0) + + for(i in 1:m) { + vids <- !is.na(xs[i,]) + nv <- sum(vids) + xv <- xs[i,vids] + + if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { + if ( !ponly ) { + vgs[i,] <- rep(NA,g) + ves[i,] <- rep(NA,g) + dfs[i,] <- rep(NA,g) + REMLs[i,] <- rep(NA,g) + stats[i,] <- rep(NA,g) + } + ps[i,] = rep(1,g) + + } + else if ( identical(x.prev, xv) ) { + if ( !ponly ) { + vgs[i,] <- vgs[i-1,] + ves[i,] <- ves[i-1,] + dfs[i,] <- dfs[i-1,] + REMLs[i,] <- REMLs[i-1,] + stats[i,] <- stats[i-1,] + } + ps[i,] <- ps[i-1,] + } + else { + if ( is.null(Z) ) { + X <- cbind(X0[vids,,drop=FALSE],xs[i,vids]) + eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X) + } + else { + vrows <- as.logical(rowSums(Z[,vids])) + X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) + eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X) + } + + for(j in 1:g) { + if ( nv == t ) { + REMLE <- emma.REMLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.R1) + if ( is.null(Z) ) { + U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),t,t,byrow=TRUE) + dfs[i,j] <- nv - q1 + } + else { + U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) + dfs[i,j] <- n - q1 + } + yt <- crossprod(U,ys[j,]) + Xt <- crossprod(U,X) + iXX <- solve(crossprod(Xt,Xt)) + beta <- iXX%*%crossprod(Xt,yt) + + if ( !ponly ) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + } + else { + if ( is.null(Z) ) { + eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids]) + nr <- sum(vids) + yv <- ys[j,vids] + REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],NULL,ngrids,llim,ulim,esp,eig.R1) + U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) + dfs[i,j] <- nr - q1 + } + else { + eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids,drop=FALSE],K[vids,vids]) + yv <- ys[j,vrows] + nr <- sum(vrows) + tv <- sum(vids) + REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],Z[vrows,vids,drop=FALSE],ngrids,llim,ulim,esp,eig.R1) + U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-tv)),nr,nr,byrow=TRUE) + dfs[i,j] <- nr - q1 + } + yt <- crossprod(U,yv) + Xt <- crossprod(U,X) + iXX <- solve(crossprod(Xt,Xt)) + beta <- iXX%*%crossprod(Xt,yt) + if (!ponly) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + } + } + ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) + } + } + } + else { + eig.L <- emma.eigen.L(Z,K) + eig.R0 <- emma.eigen.R(Z,K,X0) + + x.prev <- vector(length=0) + + for(i in 1:m) { + vids <- !is.na(xs[i,]) + nv <- sum(vids) + xv <- xs[i,vids] + + if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { + if (!ponly) { + vgs[i,] <- rep(NA,g) + ves[i,] <- rep(NA,g) + REMLs[i,] <- rep(NA,g) + dfs[i,] <- rep(NA,g) + } + ps[i,] = rep(1,g) + } + else if ( identical(x.prev, xv) ) { + if ( !ponly ) { + stats[i,] <- stats[i-1,] + vgs[i,] <- vgs[i-1,] + ves[i,] <- ves[i-1,] + REMLs[i,] <- REMLs[i-1,] + dfs[i,] <- dfs[i-1,] + } + ps[i,] = ps[i-1,] + } + else { + if ( is.null(Z) ) { + X <- cbind(X0,xs[i,]) + if ( nv == t ) { + eig.R1 = emma.eigen.R.wo.Z(K,X) + } + } + else { + vrows <- as.logical(rowSums(Z[,vids,drop=FALSE])) + X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) + if ( nv == t ) { + eig.R1 = emma.eigen.R.w.Z(Z,K,X) + } + } + + for(j in 1:g) { + vrows <- !is.na(ys[j,]) + if ( nv == t ) { + yv <- ys[j,vrows] + nr <- sum(vrows) + if ( is.null(Z) ) { + if ( nr == n ) { + REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1) + U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) + } + else { + eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE]) + REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) + } + dfs[i,j] <- nr-q1 + } + else { + if ( nr == n ) { + REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1) + U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) + } + else { + vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) + eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) + REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) + } + dfs[i,j] <- nr-q1 + } + + yt <- crossprod(U,yv) + Xt <- crossprod(U,X[vrows,,drop=FALSE]) + iXX <- solve(crossprod(Xt,Xt)) + beta <- iXX%*%crossprod(Xt,yt) + if ( !ponly ) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + } + else { + if ( is.null(Z) ) { + vtids <- vrows & vids + eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE]) + yv <- ys[j,vtids] + nr <- sum(vtids) + REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) + Xt <- crossprod(U,X[vtids,,drop=FALSE]) + dfs[i,j] <- nr-q1 + } + else { + vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids + vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE])) + eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) + yv <- ys[j,vtrows] + nr <- sum(vtrows) + REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) + U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) + Xt <- crossprod(U,X[vtrows,,drop=FALSE]) + dfs[i,j] <- nr-q1 + } + yt <- crossprod(U,yv) + iXX <- solve(crossprod(Xt,Xt)) + beta <- iXX%*%crossprod(Xt,yt) + if ( !ponly ) { + vgs[i,j] <- REMLE$vg + ves[i,j] <- REMLE$ve + REMLs[i,j] <- REMLE$REML + } + stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) + + } + } + ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) + } + } + } + if ( ponly ) { + return (ps) + } + else { + return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves)) + } +}