# R/mikecpc.R In bbolker/cpcbp: common principal components/back-projection analysis

#### Defines functions mikebppoolcpccpc_LRTcpc_objectcpcchisqcpcreorder.eehmean

```hmean <- function(L) {
## harmonic mean of a list of vectors
inv <- lapply(L,function(x) 1/x)
## compress lists into a matrix, compute means, invert ...
1/rowMeans(do.call(cbind,inv))
}

reorder.ee <- function(x) {
oo <- order(hmean(x\$evals),decreasing=TRUE)
ordfun <- function(x) {
if (is.matrix(x)) x[,oo] else x[oo]
}
## reorder all elements of all elements
## (evecs, evals) -> (by group)
## in decreasing order of harmonic mean
for (e in c("evecs","evals")) {
x[[e]] <- lapply(x[[e]], ordfun)
}
return(x)
}

cpc <- function(covs,n){
p <- nrow(covs[[1]])
ncov <- length(covs)
mlist <- list()
for(i in 1:ncov){
mlist\$evecs[[i]] <- FGalgorithm(1e-6,1e-6,p,n,covs)
}
for(j in 1:ncov){
## Flury eq. 1.20, p. 70; cov[[j]]==S_j
mlist\$evals[[j]] <- diag(t(mlist\$evecs[[j]]) %*% covs[[j]] %*% mlist\$evecs[[j]])
}
return(reorder.ee(mlist)\$evecs)
}

cpcchisq <- function(cov1,cov2,npts){
num <- sapply(cov1,det)
dem <- sapply(cov2,det)
ll <- sum((npts-1)*log(num/dem))  ## n or n-1?
return(ll)
}

cpc_object <- function(covs,npts,ncp=NULL){
if(is.null(ncp))(ncp=nrow(covs[[1]])-1)
mlist <- list()
p <- nrow(covs[[1]])
k <- length(covs)
full <- cpc(covs,npts)
if(ncp == (p-1)){
mlist\$evecs <- full
mlist\$par <- p*(p-1)/2 + k*p}
else{mlist\$evecs <- partial_cpc(covs,npts,q=ncp,B=full[[1]])
mlist\$par <- p*(p-1)/2 + k*p + (k-1)*(p-ncp)*(p-ncp-1)/2}
mlist\$evals <- list()
mlist\$cov <- list()
ncov <- length(covs)
for(j in 1:ncov){
## Flury eq. 1.20, p. 70; cov[[j]]==S_j
mlist\$evals[[j]] <- diag(t(mlist\$evecs[[j]]) %*% covs[[j]] %*% mlist\$evecs[[
j]])
mlist\$cov[[j]] <- mlist\$evecs[[j]] %*% diag(mlist\$evals[[j]])  %*%
t(mlist\$evecs[[j]])}
return(reorder.ee(mlist))
}

cpc_LRT <- function(cpc1,cpc2=NULL,covs,npts){
k <- length(covs)
p <- nrow(covs[[1]])
arbpar <- k*( p*(p-1)/2 + p)
chisqtest <- data.frame(statistic = NA,df = NA,pval = NA)
if(is.null(cpc2)){
chisqtest[1,1] <- cpcchisq(cpc1\$cov,covs,npts)
chisqtest[1,2] <- arbpar - cpc1\$par
chisqtest[1,3] <- 1 - pchisq(chisqtest[1,1],chisqtest[1,2])
}
else{
chisqtest[1,1] <- cpcchisq(cpc1\$cov,cpc2\$cov,npts)
chisqtest[1,2] <- cpc2\$par - cpc1\$par
chisqtest[1,3] <- 1 - pchisq(chisqtest[1,1],chisqtest[1,2])
}
return(chisqtest)
}

poolcpc <- function(covs,npts){
mlist <- list()
poolvar <- Reduce('+',mapply('*',covs,(npts-1),SIMPLIFY = FALSE))/(sum(npts)-length(npts))
mlist\$evals <- diag(t(eigen(poolvar)\$vectors) %*% poolvar %*% eigen(poolvar)\$vectors)
mlist\$evecs <- eigen(poolvar)\$vectors
mlist\$cov <- poolvar
return(mlist)}

mikebp <- function(x,f,evec, center = TRUE)if (missing(evec)) {
if (missing(f)) {
if (is.list(x)) {
f = x[[2]]
x = x[[1]]
}
else stop("must specify either CPC1 or a grouping factor")
}

covs <- lapply(split(data.frame(x),f),cov)
evecs = fullcpc(covs,table(f))\$evecs[[1]]
if (length(evecs) == 1 && is.na(evecs))
return(NA)
evec = evecs[, 1]
if (center)
x = scale(x, center = TRUE, scale = FALSE)
t(bpmat(evec) %*% t(x))

}
```
bbolker/cpcbp documentation built on May 10, 2017, 10:26 p.m.