library(cpca)
packageVersion("cpca")
## [1] '0.1.2'
library(plyr)
library(abind)
data(iris)
C <- daply(iris, "Species", function(x) cov(x[, -ncol(x)]))
C <- aperm(C, c(2, 3, 1)) # put the 1st dimension to the end
dim(C)
## [1] 4 4 3
dimnames(C)
## [[1]]
## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
##
## [[2]]
## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
##
## $Species
## [1] "setosa" "versicolor" "virginica"
mod <- cpc(C)
str(mod)
## List of 3
## $ D : num [1:4, 1:3] 0.1908 0.0787 0.0276 0.0121 0.4668 ...
## $ CPC : num [1:4, 1:4] 0.7467 0.4423 0.4743 0.1476 0.0911 ...
## $ ncomp: int 4
round(mod$CPC, 2)
## [,1] [,2] [,3] [,4]
## [1,] 0.75 0.09 0.63 0.20
## [2,] 0.44 -0.79 -0.33 -0.26
## [3,] 0.47 0.60 -0.54 -0.34
## [4,] 0.15 -0.02 -0.45 0.88
# See Trendafilov (2010). Stepwise estimation of common principal components.
# Computational Statistics & Data Analysis, 54(12), 3446–3457.
# doi:10.1016/j.csda.2010.03.010
# p. 10, Example 2
#
# [,1] [,2] [,3] [,4]
#[1,] 0.75 -0.09 0.63 0.20
#[2,] 0.44 0.79 -0.33 -0.26
#[3,] 0.47 -0.60 -0.54 -0.34
#[4,] 0.15 0.02 -0.45 0.88
#
# The eigenvectors must be the same, as the default method in `cpc` function
# is the power algorithm proposed by Trendafilov.
sessionInfo()
## R version 3.0.1 (2013-05-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=C LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] abind_1.4-0 plyr_1.8 cpca_0.1.2 knitr_1.3.3
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.3 evaluate_0.4.4 formatR_0.8 stringr_0.6.2
## [5] tools_3.0.1