| pca_by_als {lazy.fa} | R Documentation |
This program approximates the n x nvar input matrix Y
by the product of n x ndim F and nvar x ndim A matrices as
Y = F %*% t(A) + E
When the optimal scaling of Y is employed, the result should
be close to the one from PRINCIPALS of Young, Takane and de Leeuw (1978).
pca_by_als( Y, ndim = 2, mlevel = 4, maxiter = 1000, epsd = 1e-06, SQUAREM = 0, nSQUAREM = 1, maxalpha = -1, always = 1, reset1 = 1, reset2 = 2, print = 1, plot = 0 )
Y |
Data matrix |
ndim |
# of dimensions |
mlevel |
A vector of measurement level for each variable: |
maxiter |
Maximum # of iterations |
epsd |
The convergence criterion for the difference of param values. |
SQUAREM |
= 3 : See the help of iSQUAREM in lazy.accel package. |
nSQUAREM |
>= 1 # of iterations a which SQUAREM begins. |
maxalpha |
= -1 : See the help of iSQUAREM in lazy.accel package. |
always |
= 1 : See the help of iSQUAREM in lazy.accel package. |
reset1 |
= 0 : See the help of iSQUAREM in lazy.accel package. |
reset2 |
= 1 : See the help of iSQUAREM in lazy.accel package. |
print |
= 1 to print the result = 3 to print the intermediate result |
plot |
= 2 to plot the transformation |
nknots |
# of knots for monotone spline transformation (not available) |
minalpha |
= -999 : See the help of iSQUAREM in lazy.accel package. |
The puropose of this function is to demonstrate how the ALS (Alternating Least Squares) and inline SQUAREM can be utilized in PCA context.
The parameter set (F,A) is divided into F and A
and the conditional minimization of rss w.r.t. F and
w.r.t. A will be alternated until convergence.
When requested, the optimal scaling of each column of Y will be
performed as a part of ALS iteration.
The program may not be very efficient but it reduces rss
at each phase of iteration.
This function finds F and A matrices of rank ndim
and the optimal transformation of each column of the data matrix Y
which minimize
sum over j of( (Yhat[,j]-(F%*%t(A))[,j])^2 / var(Yhat[,j]) )
where Yhat=t(Y) is the column-wise optimally scaled Y.
The actual minimization is done by normalizing Yhat so that
var(Yhat[,j])=1 and minimizing
rss=sum( (Yhat-F%*%A)^2 )
by alternating least squares..
The core part of this function is as follows:
# Repeat the following 4 phases until convergence. # 1 Given Yhat and F, estimate A A <- t( solve(t(F)%*%F)%*%t(F)%*%Yhat ) # 2 Given Yhat and F, estimate F F <- t( solve(t(A)%*%A)%*%t(A)%*%t(Yhat) ) # 3 Given F and A, estimate each column of Yhat Yhat[,j] <- optimally scaled Y[,j] which minimizes rss. # 4 Column standardize Yhat. # Check convergence: rss should decrease during the iteration. #
Note that the column standardized Y will be used as the initial
Yhat.
Note that when mlevel=0 rss will be minimized without
normalization of Yhat. Therefore, the result should be equievalent
to PCA after suitable normalization of F and A.
A list of
F n x ndim componet score matrix whose dispersion matrix is I A nvar x ndim componet loadings matrix satisfying t(A)%*%A = diagonal Yhat n x nvar optimally scaled data matrix rmse = sqrt( rss/(n*nvar) ) mlevel transformation level
Young, W. F., Takane, Y. and de Leeuw, J. (1978)
The principal components of mixed measurement level multivariate data:
An alternating least squares method with optimal scaling features.
Psychometrika volume 43, pages 279-281
# # Comparison of pca by als and native prcomp # # # generate column centered Y seed <- 1701 set.seed(seed) n <- 50; nvar <- 6 Y <- matrix(rnorm(n*nvar),n,nvar) Y <- scale( Y, center=TRUE ) attributes(Y)$`scale.center`=NULL # pca by als (no optimal transformation) ndim <- 2 res <- pca_by_als( Y, ndim, SQUAREM=0, always=1, print=1, mlevel=0 ) # normalize the result so that Disp(F)=I and t(A)%*%A=diagonal FA <- normalize_fa( res$F, res$A ) F1 <- FA$F; A1 <- FA$A rmse1 <- sqrt(ssq(Y-F1%*%t(A1))/n/nvar) # pca by native prcomp: Y == F2 %*% t(A2), t(F2)%*%F2 = diag, t(A2)%*%A2 = I resprc <- prcomp( Y, rank=ndim ) sqrteva <- resprc$sdev[1:ndim]*sqrt((n-1)) # vardef adjusted, eva of Y'Y A2 <- resprc$rotation # get F2 by solving Y=F2%*%t(A2) or t(Y)=A2%*%t(F2) for F2 F2 <- Y%*%A2%*%solve(t(A2)%*%A2) # from prcomp to pca: # Note that sqrt eva from native prcomp is the eva of Y'Y, not Y'Y/n A21 <- A2%*%diag(sqrteva)/sqrt(n) F21 <- F2%*%diag(sqrt(n)/sqrteva) rmse21 <- sqrt(ssq(Y-F21%*%t(A21))/n/nvar) # FA <- normalize_fa( F21, A21 ) # F21n <- FA$F; A21n <- FA$A Print(A1,A21) Print(rmse1,rmse21) # Simulating PRINCIPALS # generate F according to the FA model seed <- 1701 set.seed(seed) n <- 200; nvar <- 9 ndim0 <- 3 F <- matrix(rnorm(n*ndim0),n,ndim0) A <- matrix(rnorm(nvar*ndim0), nvar,ndim0) A[,3] <- 0.5*A[,3] # transform Y Y <- F%*%t(A) + 0.2*matrix(rnorm(n*nvar),n,nvar) # transform Y Y[,2]=(Y[,2]-min(Y[,2]))^3 Y[,3]=(Y[,3]-min(Y[,3]))^3 Y[,4]=(Y[,4]-min(Y[,4]))^(1/3) Y[,5]=(Y[,5]-min(Y[,5]))^(1/3) Y <- scale(Y, center=TRUE, scale=TRUE) # round Y Y <- round(10*Y,0) # PRINCIPALS ndim <- 2 res5 <- pca_by_als(Y, ndim, SQUAREM=3, always=1, print=1, plot=2, mlevel=5) res4 <- pca_by_als(Y, ndim, SQUAREM=3, always=1, print=1, plot=2, mlevel=4) res3 <- pca_by_als(Y, ndim, SQUAREM=3, always=1, print=1, plot=2, mlevel=3) Print(res5$mlevel, res5$rmse) Print(res4$mlevel, res4$rmse) Print(res3$mlevel, res3$rmse)