R/Mpca.R

Mpca <-
function(x,N1=3,N2=2,tol=.001,N2p=10,Nran=50,
Nkeep=10,SEED=TRUE,op.pro=.1,SCORES=FALSE,pval=NULL){
#
# Robust PCA using Marrona's method (2005, Technometrics)
#
# x is an N by m matrix containing data
# N1, N2, N2p, Nran and Nkeep indicate how many
# iterations are used in the various portions of the
# Marrona robust PCA; see Marrona's paper.
#
# op.pro is the maximum  proportion of unexplained
#   variance that is desired. If pval is not specified, will
#  add variables until this proportion is less than op.pro.
#
# pval, if specified, will use p=pval of the m variables only and report
#  the proportion of unexplained variance.
#  The weighted covariance matrix is returned as well.
#
# SCORES=T, scores are reported and return based on the number of
# variables indicated by pval. pval must be specified.
#
# pval not specified, computes proportion of unexplained  variance
# using p=1, 2 ... variables; results returned in
#
scores<-NULL
wt.cov<-NULL
x<-elimna(x)
if(SEED)set.seed(2)
m<-ncol(x)
n<-nrow(x)
bot<-marpca(x,p=0,N1=N1,N2=N2,tol=tol,N2p=N2p,Nran=Nran,Nkeep=Nkeep,SEED=SEED)
bot<-bot$var.op
mn1<-m-1
rat<-1
it<-0
ratval<-NULL
if(is.null(pval)){
ratval<-matrix(nrow=mn1,ncol=2)
dimnames(ratval)<-list(NULL,c("p","pro.unex.var"))
ratval[,1]<-c(1:mn1)
for(it in 1:mn1){
if(rat>op.pro){
temp<-marpca(x,p=it,N1=N1,N2=N2,tol=tol,N2p=N2p,Nran=Nran,Nkeep=Nkeep,
SEED=SEED)
rat<-temp$var.op/bot
ratval[it,2]<-rat
}}}
if(!is.null(pval)){
if(pval>=m)stop("This method assumes pval<ncol(x)")
temp<-marpca(x,p=pval,N1=N1,N2=N2,tol=tol,N2p=N2p,Nran=Nran,Nkeep=Nkeep,
SEED=SEED)
wt.cov<-temp$wt.cov
ratval<-temp$var.op/bot
}
if(SCORES){
if(is.null(pval))stop("When computing scores, need to specify pval")
temp2<-marpca(x,ncol(x))
ev<-eigen(temp2$wt.cov)
ord.val<-order(ev$values)
mn1<-m-pval+1
Bp<-ev$vectors[,ord.val[mn1:m]] #m by m
xmmu<-x
for(j in 1:m)xmmu[,j]<-x[,j]-temp2$wt.mu[j]
scores<-matrix(ncol=pval,nrow=n)
for(i in 1:n)scores[i,]<-t(Bp)%*%as.matrix(xmmu[i,])
}
list(B=temp$B,a=temp$a,var.op=temp$var.op,unexplained.pro.var=ratval,
scores=scores,wt.cov=wt.cov)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.