# Purpose: Repeated Score test for multivariate outcome regression
# Updated: 180910
########################
# Multivariate case
########################
#' Multivariate Regression Repeated Score Test
#'
#' Individually tests each column of \code{G} for association with the target
#' outcome. \code{G} may contain missing elements, although the remaining model
#' matrices may not. \code{rScore.mnr} accelerates association testing by
#' recycling the same null model for each hypothesis test.
#'
#' @param Y Outcome matrix.
#' @param j Column number of the outcome of interest. By default, \code{j=1}.
#' @param G Numeric matrix of covariates for the target outcome whose regression
#' coefficients are zero under the null.
#' @param X List of model matrices, one for each outcome.
#' @param maxit Maximum number of parameter updates.
#' @param eps Minimum acceptable improvement in log likelihood.
#' @param report Report model fitting progress? Default is FALSE.
#' @param parallel Run association testing in parallel? Must register parallel
#' backend first.
#'
#' @importFrom foreach foreach '%do%'
#' @importFrom plyr aaply
#' @importFrom stats model.matrix pchisq resid vcov
#' @export
#' @return A numeric vector of p-values, one for each column of \code{G}.
#' @examples
#' \dontrun{
#' # See `? rMNR` for data generation
#' # See vignette for test description
#' # Genotype matrix
#' G = replicate(2000,rbinom(n=1000,size=2,prob=0.25));
#' storage.mode(G) = "numeric";
#' # Introduce missingness
#' G[sort(sample(length(G),size=0.01*length(G)))] = NA;
#' # Repeated Score test
#' R = rScore.mnr(Y=Y,G=G,X=X,report=T);
#' # Estimated size
#' mean(R<=0.05);
#' }
rScore.mnr = function(Y,j=1,G,X,maxit=100,eps=1e-8,report=F,parallel=F){
# Check input type
if(!is.matrix(Y)){stop("A numeric matrix is expected for Y.")};
if(!is.matrix(G)){stop("A numeric matrix is expected for G.")}
if(!is.list(X)||(!all(unlist(lapply(X,is.matrix))))){
stop("A list of numeric model matrices is expected for X.")};
# Restructure
k = ncol(Y);
# Place target outcome first
perm = c(j,seq(from=1,to=k)[-j]);
Y = Y[,perm];
X = X[perm];
rm(j);
# Check for missingness
Miss = sum(is.na(Y))+sum(unlist(lapply(X,function(x){sum(is.na(x))})));
if(Miss>0){stop("Inputs other than G should contain no missing data.")};
# Fit null model
M0 = fit.mnr(Y=Y,X=X,maxit=maxit,eps=eps,report=report);
# Extract covariance
L = vcov(M0,type="Outcome",inv=T);
# Extract residuals
E = resid(M0);
# Working vector
u = MMP(E,L[1,]);
# Function to calculate score statistics
aux = function(g){
# Check for missingness
key = !is.na(g);
# If missingness, exclude
if(sum(!key)>0){
u0 = u[key];
g0 = g[key];
X0 = lapply(X,function(x){x[key,,drop=F]});
X1 = lapply(X,function(x){x[!key,,drop=F]});
} else {
u0 = u;
g0 = g;
X0 = X;
}
## Information matrices
# Marginal information
Ibb = L[1,1]*matIP(g0,g0);
i = NULL;
Iba = foreach(i=1:k,.combine=cbind) %do% {
return(matIP(g0,X0[[i]])*L[1,i])
};
# Nusiance information
Iaa = vcov(M0,type="Information");
if(sum(!key)>0){
l = NULL;
# Loss of information
Jaa = foreach(i=1:k,.combine=rbind) %do% {
Slice = foreach(l=1:k,.combine=cbind) %do% {
Out = matIP(X1[[i]],X1[[l]])*L[i,l];
return(Out);
};
};
Iaa = Iaa-Jaa;
};
# Efficient information
V = as.numeric(SchurC(Ibb=Ibb,Iaa=Iaa,Iba=Iba));
# Score vector
U = as.numeric(matIP(g0,u0));
# Score statistic
Ts = (U^2)/V;
return(Ts);
}
# Score statistics
S = aaply(.data=G,.margins=2,.fun=aux,.parallel=parallel);
# P values
p = pchisq(q=S,df=1,lower.tail=F);
return(p);
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.