lpls <-
function(X1,X2,X3,npc=2,doublecenter=TRUE,scaledata=c(FALSE,FALSE,FALSE),
type=c("exo"), impute=FALSE, niter=25, subsetX1=NULL, subsetX3=NULL,...){
#The three matrices are assumed positioned in this manner
# _________
# | |
# | |
# | t(X3) |
# | |
# |_________|
#
#
# _______ _________
# | | | |
# | | | |
# | X1 | | X2 |
# | | | |
# |_______| |_________|
#
# if(type!="exo_ort") impute=TRUE
if(impute){
if(any(is.na(X1))) X1 <-.svd.imp(X1)
if(any(is.na(X2))) X2 <-.svd.imp(X2)
if(any(is.na(X3))) X3 <-.svd.imp(X3)
}
#Extracting subsets
if(!is.null(subsetX1)){
X1 <- X1[subsetX1,,drop=F]
X2 <- X2[subsetX1,,drop=F]
}
if(!is.null(subsetX3)){
X3 <- X3[subsetX3,,drop=F]
X2 <- X2[,subsetX3,drop=F]
}
#For saving
X1save <- X1
X2save <- X2
X3save <- X3
#Dimensions
X1dim<-dim(X1)
X2dim<-dim(X2)
X3dim<-dim(X3)
#Centering and scaling
X1 <- scale(X1,scale=scaledata[1])
mX1 <- attr(X1,"scaled:center")
X3 <- scale(X3,scale=scaledata[3])
mX3 <- attr(X3,"scaled:center")
rowmX2 <- apply(X2,1,mean, na.rm=TRUE)
colmX2 <- apply(X2,2,mean, na.rm=TRUE)
grandmX2 <- mean(X2, na.rm=TRUE)
#Do double centering of X2?
if(doublecenter){
X2 <- X2 - t(matrix(rep(1,X2dim[2]),ncol=1)%*%rowmX2) -
matrix(rep(1,X2dim[1]),ncol=1)%*%colmX2 +
matrix(grandmX2,nrow=X2dim[1],ncol=X2dim[2])
}else{
X2 <- X2 - matrix(grandmX2,nrow=X2dim[1],ncol=X2dim[2])
}
#Column scaling of X2
X2 <- scale(X2,scale=scaledata[2])
attributes(X1save) <- attributes(X1)
attributes(X2save) <- attributes(X2)
attributes(X3save) <- attributes(X3)
dlist <- list(X1=as.matrix(t(X1)), X2=as.matrix(X2), X3=as.matrix(X3))
T11 <- T12 <- T21 <- T22 <- T31 <- T32 <- numeric(0)
W21 <- W22 <- numeric(0)
P1 <- P3 <- P21 <- P22 <- numeric(0)
D <- diag(rep(1,npc))
X1totvar <- X2totvar <- X3totvar <- rep(0,npc+1)
for(a in 1:npc){
X1totvar[a] <-sum(dlist$X1^2, na.rm=TRUE)
X2totvar[a] <-sum(dlist$X2^2, na.rm=TRUE)
X3totvar[a] <-sum(dlist$X3^2, na.rm=TRUE)
latent <- .extractscores(dlist, niter=niter)
#Score vectors
T11<-cbind(T11, .lnorm(latent[[1]]$t1))
T12<-cbind(T12, .lnorm(latent[[1]]$t2))
T21<-cbind(T21, .lnorm(latent[[2]]$t1))
T22<-cbind(T22, .lnorm(latent[[2]]$t2))
T31<-cbind(T31, .lnorm(latent[[3]]$t1))
T32<-cbind(T32, .lnorm(latent[[3]]$t2))
#Fitting an exo-LPLS model
if(type=="exo_ort"){
#Orthogonal scores
P1 <- cbind(P1, .projectonto(dlist$X1,T21[,a]))
P3 <- cbind(P3, .projectonto(dlist$X3,T22[,a]))
P21 <- cbind(P21, .projectonto(dlist$X2,T21[,a]))
P22 <- cbind(P22, .projectonto(dlist$X2,T22[,a]))
d <- drop(solve(crossprod(T21[,a]))%*%t(T21[,a])%*%P22[,a])
D[a,a] <- d
#Deflation
dlist$X1 <- dlist$X1 - P1[,a]%*%t(T21[,a])
dlist$X3 <- dlist$X3 - T22[,a]%*%t(P3[,a])
dlist$X2 <- dlist$X2 - T21[,a]%*%t(P21[,a]) - P22[,a]%*%t(T22[,a]) + T21[,a]%*%t(T22[,a])*drop(d)
}#end type=="exo"
if(type=="exo"){
#Simpler model for X2, but non-orthogonal scores
P1 <- .projectonto(t(X1),T21)
P3 <- .projectonto(t(X3),T22)
P21 <- .projectonto(t(X2),T21)
P22 <- .projectonto(X2,T22)
D <- solve(crossprod(T21))%*%t(T21)%*%P22
#Deflation
dlist$X1 <- t(X1 - T21%*%t(P1))
dlist$X3 <- X3 - T22%*%t(P3)
dlist$X2 <- X2 - T21%*%D%*%t(T22)
}#end type=="exo"
if(type=="endo"){
P1 <- .projectonto(t(X1),T12)
P3 <- .projectonto(t(X3),T31)
D <- solve(crossprod(T12))%*%(t(T12)%*%.projectonto(X2, T31))
#Deflation
dlist$X1 <- t(X1 - T12%*%t(P1))
dlist$X3 <- X3 - T31%*%t(P3)
dlist$X2 <- X2 - T12%*%(D%*%t(T31))
}#end type=="endo"
}
X1totvar[npc+1] <-sum(dlist$X1^2, na.rm=TRUE)
X2totvar[npc+1] <-sum(dlist$X2^2, na.rm=TRUE)
X3totvar[npc+1] <-sum(dlist$X3^2, na.rm=TRUE)
X1varprop <- diff((X1totvar[1]-X1totvar)/(X1totvar[1]))
X2varprop <- diff((X2totvar[1]-X2totvar)/(X2totvar[1]))
X3varprop <- diff((X3totvar[1]-X3totvar)/(X3totvar[1]))
#Various output
B1 <- B3 <- Ca <- NULL
if(type!="endo"){
B1 <- T31%*%(solve(crossprod(P21,T31))%*%t(P1))
B3 <- T12%*%(solve(crossprod(P22,T12))%*%t(P3))
options(warn=-1)
R1 <- cor(X1,T21, use="pairwise.complete.obs")
R21 <- t(cor(T21,X2, use="pairwise.complete.obs"))
R22 <- t(cor(T22,t(X2), use="pairwise.complete.obs"))
R3 <- cor(X3,T22, use="pairwise.complete.obs")
R2rmean <- cor(rowmX2,T21, use="pairwise.complete.obs")
R2cmean <- cor(colmX2,T22, use="pairwise.complete.obs")
options(warn=0)
}else if(type=="endo"){
V1 <- T11%*%(solve(crossprod(P1,T11)))
V3 <- T32%*%(solve(crossprod(P3,T32)))
Ca <- V1%*%D%*%t(V3)
options(warn=-1)
R1 <- cor(X1,T12, use="pairwise.complete.obs")
R21 <- t(cor(T12,X2, use="pairwise.complete.obs"))
R22 <- t(cor(T31,t(X2), use="pairwise.complete.obs"))
R3 <- cor(X3,T31, use="pairwise.complete.obs")
R2rmean <- cor(rowmX2,T12, use="pairwise.complete.obs")
R2cmean <- cor(colmX2,T31, use="pairwise.complete.obs")
options(warn=0)
}
res <- list(call=match.call())
res$ncomp <- npc
res$coefficients <- list(B1=B1, B3=B3, C=Ca)
res$scores <- list(T11=T11, T12=T12, T21=T21, T22=T22, T31=T31, T32=T32)
res$loadings <- list(P1=P1, P3=P3, P21=P21, P22=P22, D=D)
res$corloadings <- list(R1=R1, R21=R21, R22=R22, R3=R3, R2rmean=R2rmean, R2cmean=R2cmean)
res$means <- list(mX1=mX1, mX3=mX3, grandmX2=grandmX2, rowmX2=rowmX2, colmX2=colmX2)
res$data <- list(X1=X1save, X2=X2save, X3=X3save)
res$residuals <- dlist
res$options <- list(doublecenter=doublecenter, scaledata=scaledata, type=type)
res$vars <- list(X1varprop=X1varprop, X2varprop=X2varprop, X3varprop=X3varprop)
class(res)<-c("lpls")
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.