# R/matrices.R In kkholst/lava: Latent Variable Models

#### Defines functions matrices.lvmmatrices2matrices.multigroupmatrices.lvmmat.lvm`matrices`

```###{{{

`matrices` <-
function(x,...) UseMethod("matrices")

###{{{ matrices.lvm

mat.lvm <- function(x,ii=index(x),...) {
A <- ii\$A ## Matrix with fixed parameters and ones where parameters are free
M1 <- ii\$M1 ## Index of free and _unique_ regression parameters
P <- ii\$P  ## Matrix with fixed variance parameters and ones where parameters are free
P1 <- ii\$P1 ## Index of free and _unique_ regression parameters

constrain.par <- names(constrain(x))
parval <- list()

parBelongsTo <- list(mean=seq_len(ii\$npar.mean),
reg=seq_len(ii\$npar.reg)+ii\$npar.mean,
cov=seq_len(ii\$npar.var)+ii\$npar.mean+ii\$npar.reg,
epar=seq_len(ii\$npar.ex)+with(ii,npar.reg+npar.var+npar.mean),
cpar=numeric())

idxA <- which(M1==1)
pidxA <- parBelongsTo\$reg
if (ii\$npar.reg>0) {
A[idxA] <- pidxA
for (p in ii\$parname) {
idx <- which((x\$par==p))
newval <- A[idx[1]]
attributes(newval)\$reg.idx <- idx
attributes(newval)\$reg.tidx <- which(t(x\$par==p))
parval[[p]] <- newval
if (length(idx)>1) {
idxA <- c(idxA,idx[-1])
pidxA <- c(pidxA,rep(A[idx[1]],length(idx)-1))
A[idx] <- A[idx[1]]
}
} ## duplicate parameters
}

pars.var <- parBelongsTo\$cov
idxdiag <- (seq(ncol(P1))-1)*ncol(P1) + seq(ncol(P1))
idxP <- idxdiag[which(P1[idxdiag]==1)]
pidxP <- pars.var[seq_len(length(idxP))]
P[idxP] <- pidxP

pars.off.diag <- pars.var
if (length(pidxP)>0) {
pars.off.diag <- pars.off.diag[-seq_len(length(pidxP))]
}

counter <- 0
if (length(pars.off.diag)>0 & ncol(P)>1)
for (i in seq_len(ncol(P1)-1))
for (j in seq(i+1,nrow(P1))) {
if (ii\$P1[j,i]!=0) {
counter <- counter+1
pos <- c(j+(i-1)*ncol(P1),
i+(j-1)*ncol(P1))
P[j,i] <- P[i,j] <- pars.off.diag[counter]
idxP <- c(idxP,pos); pidxP <- c(pidxP,P[j,i],P[i,j])
}
}

if (length(ii\$covparname)>0)
for (p in ii\$covparname) {
idx <- which(x\$covpar==p)
isOffDiag <- !(idx[1]%in%idxdiag)
if (!(p%in%ii\$parname)) {
parval[[p]] <- P[idx[1]]
}
attributes(parval[[p]])\$cov.idx <- idx
if (length(idx)>1+isOffDiag) {
P[idx[-seq(1+isOffDiag)]] <- parval[[p]]
}
if (ii\$npar.reg>0 && p%in%ii\$parname) {
parBelongsTo\$reg <- c(parBelongsTo\$reg,p)
idx.reg <- which(x\$par==p)
P[idx] <- A[idx.reg[1]]
atr <- attributes(parval[[p]])
parval[[p]] <- A[idx.reg[1]]
attributes(parval[[p]]) <- atr
idxP <- c(idxP,idx)
pidxP <- c(pidxP,rep(P[idx[1]],length(idx)))
} else {
idxP <- c(idxP,idx[-seq(1+isOffDiag)])
pidxP <- c(pidxP,rep(P[idx[1]],length(idx)-1-isOffDiag))
}
} ## duplicate parameters

idxM <- c()
pidxM <- seq_len(ii\$npar.mean)
v <- NULL

## named <- sapply(x\$mean, function(y) is.character(y) & !is.na(y))
fixed <- sapply(x\$mean, function(y) is.numeric(y) & !is.na(y))
v <- rep(0,length(x\$mean))
names(v) <- colnames(P)
if (ii\$npar.mean>0) {
idxM <- which(ii\$v1==1)
v[idxM] <- pidxM
}
if (any(fixed))
v[fixed] <- unlist(x\$mean[fixed])

for (p in ii\$mparname) {
idx <- which(x\$mean==p)
if (!(p%in%c(ii\$parname,ii\$covparname))) {
if (length(idx)>1) {
pidxM <- c(pidxM,rep(v[idx[1]],length(idx)-1))
idxM <- c(idxM,idx[-1])
}
parval[[p]] <- v[idx[1]]
v[idx] <- parval[[p]]
}
attributes(parval[[p]])\$m.idx <- idx
if (p %in% ii\$covparname & !(p %in% ii\$parname)) {
parBelongsTo\$cov <- c(parBelongsTo\$cov,p)
idx.2 <- which(x\$covpar==p)
v[idx] <- P[idx.2[1]]
pidxM <- c(pidxM,rep(P[idx.2[1]],length(idx)))
idxM <- c(idxM,idx)
}
if (p %in% ii\$parname) {
parBelongsTo\$reg <- c(parBelongsTo\$reg,p)
idx.2 <- which(x\$par==p)
v[idx] <- A[idx.2[1]]
pidxM <- c(pidxM,rep(A[idx.2[1]],length(idx)))
idxM <- c(idxM,idx)
}
}

## Ex-parameters
idxE <- NULL
pidxE <- parBelongsTo\$epar
## named <- sapply(x\$exfix, function(y) is.character(y) & !is.na(y))
fixed <- sapply(x\$exfix, function(y) is.numeric(y) & !is.na(y))
epar <- rep(0,length(x\$exfix))
names(epar) <- names(x\$expar)
if (!(ii\$npar.ex==0)) {
idxE <- which(ii\$e1==1)
epar[idxE] <- pidxE
}
if (any(fixed))
epar[fixed] <- unlist(x\$exfix[fixed])
for (p in ii\$eparname) {
idx <- which(x\$exfix==p)
if (!(p%in%c(ii\$parname,ii\$covparname,ii\$mparname))) {
if (length(idx)>1) {
idxE <- c(idxE,idx[-1])
pidxE <- c(pidxE,rep(epar[idx[1]],length(idx)-1))
}
parval[[p]] <- epar[idx[1]]
}
attributes(parval[[p]])\$e.idx <- idx

if (length(idx)>1)
epar[idx[-1]] <- parval[[p]]
if (p %in% setdiff(ii\$covparname,c(ii\$parname,ii\$mparname))) {
parBelongsTo\$cov <- c(parBelongsTo\$cov,p)
idx.2 <- which(x\$covpar==p)
epar[idx] <- P[idx.2[1]]
pidxE <- c(pidxE,rep(P[idx.2[1]],length(idx)))
idxE <- c(idxE,idx)
}
if (p %in% setdiff(ii\$parname,ii\$mparname)) {
parBelongsTo\$reg <- c(parBelongsTo\$reg,p)
idx.2 <- which(x\$par==p)
epar[idx] <- A[idx.2[1]]
pidxE <- c(pidxE,rep(A[idx.2[1]],length(idx)))
idxE <- c(idxE,idx)
}
if (p %in% ii\$mparname) {
parBelongsTo\$mean <- c(parBelongsTo\$mean,p)
idx.2 <- which(x\$mean==p)
epar[idx] <- v[idx.2[1]]
pidxE <- c(pidxE,rep(v[idx.2[1]],length(idx)))
idxE <- c(idxE,idx)
}
}
ee <- cbind(idxE,pidxE); rownames(ee) <- names(x\$expar)[ee[,1]]

## Constrained...
constrain.par <- names(constrain(x))
constrain.idx <- NULL

if (length(constrain.par)>0) {
constrain.idx <- list()
for (p in constrain.par) {
reg.tidx <- reg.idx <- cov.idx <- m.idx <- e.idx <- NULL
myc <- constrain(x)[[p]]
xargs <- manifest(x)[na.omit(match(attributes(myc)\$args,manifest(x)))]
if (length(xargs)>0) {
parval[xargs] <- 0
}
if (p%in%ii\$parname.all) {
reg.idx <- which(x\$par==p)
reg.tidx <- which(t(x\$par==p))
}
if (p%in%ii\$covparname.all) {
cov.idx <- which(x\$covpar==p)
}
if (p%in%ii\$mparname.all) {
m.idx <- which(x\$mean==p)
}
if (p%in%ii\$eparname.all) {
e.idx <- which(x\$exfix==p)
}
constrain.idx[[p]] <- list(reg.idx=reg.idx,reg.tidx=reg.tidx,cov.idx=cov.idx,m.idx=m.idx,e.idx=e.idx)
}
}

parBelongsTo <- lapply(parBelongsTo,function(x) sort(unique(x)))

return(list(mean=cbind(idxM,pidxM),
reg=cbind(idxA,pidxA),
cov=cbind(idxP,pidxP),
epar=ee,
parval=parval,
constrain.idx=constrain.idx,
parBelongsTo=parBelongsTo))

}

matrices.lvm <- function(x,pars,meanpar=NULL,epars=NULL,data=NULL,...) {
ii <- index(x)
pp <- c(rep(NA,ii\$npar.mean),pars,epars)
##v <- NULL
v <- ii\$v0
if (!is.null(meanpar) && length(meanpar)>0) {
pp[seq(ii\$npar.mean)] <- meanpar
v[ii\$mean[,1]] <- meanpar[ii\$mean[,2]]
}

A <- ii\$A
A[ii\$reg[,1]] <- pp[ii\$reg[,2]]
P <- ii\$P
P[ii\$cov[,1]] <- pp[ii\$cov[,2]]
e <- NULL

if (length(x\$expar)>0) {
e <- rep(0,length(x\$expar))
fixed <- sapply(x\$exfix, function(y) is.numeric(y) & !is.na(y))
if (any(fixed))
e[fixed] <- unlist(x\$exfix[fixed])
if (nrow(ii\$epar)>0)
e[ii\$epar[,1]] <- pp[ii\$epar[,2]]
names(e) <- names(x\$expar)
}

parval <- lapply(ii\$parval,function(x) {
res <- pp[x];
attributes(res) <- attributes(x);
res })
## Constrained...
constrain.par <- names(constrain(x))
constrain.idx <- NULL
cname <- constrainpar <- c()
if (length(constrain.par)>0 && is.numeric(c(pars,meanpar))) {
constrain.idx <- list()
for (p in constrain.par) {
cname <- c(cname,p)
myc <- constrain(x)[[p]]
xargs <- manifest(x)[na.omit(match(attributes(myc)\$args,manifest(x)))]
if (length(xargs)>0) {
if (!is.null(data)) {
parval[xargs] <- (data)[xargs]
} else parval[xargs] <- 0
}
val <- unlist(c(parval,constrainpar,x\$mean,e)[attributes(myc)\$args])
cpar <- myc(val);
constrainpar <- c(constrainpar,list(cpar)); names(constrainpar) <- cname
if (p%in%ii\$parname.all) {
if (!is.null(val))
A[ii\$constrain.idx[[p]]\$reg.idx] <- cpar
}
if (p%in%ii\$covparname.all) {
if (!is.null(val))
P[ii\$constrain.idx[[p]]\$cov.idx] <- cpar
}
if (p%in%ii\$mparname.all) {
if (!is.null(val))
v[ii\$constrain.idx[[p]]\$m.idx] <- cpar
}
if (p%in%ii\$eparname.all) {
if (!is.null(val))
e[ii\$constrain.idx[[p]]\$e.idx] <- cpar
}
}
}

return(list(A=A, P=P, v=v, e=e, parval=parval,
constrain.idx=ii\$constrain.idx, constrainpar=constrainpar))
}

###}}} matrices.lvm

###{{{ matrices.multigroup

matrices.multigroup <- function(x, p, ...) {
pp <- modelPar(x,p)
res <- list()
for (i in seq_len(x\$ngroup))
res <- c(res, list(matrices2(x\$lvm[[i]],pp\$p[[i]])))
return(res)
}

###}}}

matrices2 <- function(x,p,...) {
m0 <- p[seq_len(index(x)\$npar.mean)]
p0 <- p[with(index(x),seq_len(npar)+npar.mean)]
e0 <- p[with(index(x),seq_len(npar.ex)+npar.mean+npar)]
matrices(x,p0,m0,e0,...)
}

###{{{ matrices, to be superseeded by above definition

matrices.lvm <- function(x,pars,meanpar=NULL,epars=NULL,data=NULL,...) {
ii <- index(x)
A <- ii\$A ## Matrix with fixed parameters and ones where parameters are free
M1 <- ii\$M1 ## Index of free and _unique_ regression parameters
P <- ii\$P  ## Matrix with fixed variance parameters and ones where parameters are free
P1 <- ii\$P1 ## Index of free and _unique_ regression parameters

constrain.par <- names(constrain(x))
parval <- list()

if (ii\$npar.reg>0) {
A[which(M1==1)] <- pars[seq_len(ii\$npar.reg)]
for (p in ii\$parname) {
idx <- which((x\$par==p))
newval <- A[idx[1]]
attributes(newval)\$reg.idx <- idx
attributes(newval)\$reg.tidx <- which(t(x\$par==p))
parval[[p]] <- newval
if (length(idx)>1) {
A[idx[-1]] <- parval[[p]]
}
} ## duplicate parameters
}

if (ii\$npar.reg==0) {
pars.var <- pars
} else {
pars.var <- pars[-seq_len(ii\$npar.reg)]
}

diag(P)[ii\$which.diag] <- pars.var[seq_along(ii\$which.diag)]
pars.off.diag <- pars.var
if (length(ii\$which.diag)>0)
pars.off.diag <- pars.off.diag[-seq_along(ii\$which.diag)]
counter <- 0
if (length(pars.off.diag)>0 & ncol(P)>1)
for (i in seq_len(ncol(P1)-1))
for (j in seq(i+1,nrow(P1))) {
if (ii\$P1[j,i]!=0) {
counter <- counter+1
P[j,i] <- pars.off.diag[counter]
}
}

if (length(ii\$covparname)>0)
for (p in ii\$covparname) {
idx <- which(x\$covpar==p)
if (!(p%in%ii\$parname)) {
parval[[p]] <- P[idx[1]]
}
attributes(parval[[p]])\$cov.idx <- idx
if (length(idx)>1) {
P[idx[-1]] <- parval[[p]]
}
if (ii\$npar.reg>0 && p%in%ii\$parname) {
idx.reg <- which(x\$par==p)
P[idx] <- A[idx.reg[1]]
atr <- attributes(parval[[p]])
parval[[p]] <- A[idx.reg[1]] ###?????
attributes(parval[[p]]) <- atr
}
} ## duplicate parameters
P[upper.tri(P)] <- t(P)[upper.tri(P)]  ## Symmetrize...

v <- NULL
{
## named <- sapply(x\$mean, function(y) is.character(y) & !is.na(y))
fixed <- sapply(x\$mean, function(y) is.numeric(y) & !is.na(y))
v <- rep(0,length(x\$mean))
names(v) <- colnames(P)
if (!(is.null(meanpar) | ii\$npar.mean==0))
v[ii\$v1==1] <- meanpar
if (any(fixed))
v[fixed] <- unlist(x\$mean[fixed])

for (p in ii\$mparname) {
idx <- which(x\$mean==p)

if (!(p%in%c(ii\$parname,ii\$covparname))) {
parval[[p]] <- v[idx[1]]
}
attributes(parval[[p]])\$m.idx <- idx

if (length(idx)>1)
v[idx[-1]] <- parval[[p]]
if (p %in% ii\$covparname & !(p %in% ii\$parname)) {
idx.2 <- which(x\$covpar==p)
v[idx] <- P[idx.2[1]]
}
if (p %in% ii\$parname) {
idx.2 <- which(x\$par==p)
v[idx] <- A[idx.2[1]]
}
}
}

## Ex-parameters
e <- NULL
{
## named <- sapply(x\$exfix, function(y) is.character(y) & !is.na(y))
fixed <- sapply(x\$exfix, function(y) is.numeric(y) & !is.na(y))

e <- rep(0,length(x\$exfix))

names(e) <- names(x\$expar)
if (!(is.null(epars) | ii\$npar.ex==0))
e[which(ii\$e1==1)] <- epars
if (any(fixed))
e[fixed] <- unlist(x\$exfix[fixed])
for (p in ii\$eparname) {
idx <- which(x\$exfix==p)
if (!(p%in%c(ii\$parname,ii\$covparname,ii\$mparname))) {
parval[[p]] <- e[idx[1]]
}
attributes(parval[[p]])\$e.idx <- idx

if (length(idx)>1)
e[idx[-1]] <- parval[[p]]
if (p %in% setdiff(ii\$covparname,c(ii\$parname,ii\$mparname))) {
idx.2 <- which(x\$covpar==p)
e[idx] <- P[idx.2[1]]
}
if (p %in% setdiff(ii\$parname,ii\$mparname)) {
idx.2 <- which(x\$par==p)
e[idx] <- A[idx.2[1]]
}
if (p %in% ii\$mparname) {
idx.2 <- which(x\$mean==p)
e[idx] <- v[idx.2[1]]
}
}
}

## Constrained...
constrain.idx <- NULL
cname <- constrainpar <- c()
if (length(constrain.par)>0 && is.numeric(c(pars,meanpar,e))) {
constrain.idx <- list()
for (p in constrain.par) {
cname <- c(cname,p)
reg.tidx <- reg.idx <- cov.idx <- m.idx <- e.idx <- NULL
myc <- constrain(x)[[p]]
xargs <- manifest(x)[na.omit(match(attributes(myc)\$args,manifest(x)))]
if (length(xargs)>0) {
if (!is.null(data)) {
parval[xargs] <- (data)[xargs]
} else parval[xargs] <- 0
}
val <- rbind(unlist(c(parval,constrainpar,x\$mean,e)[attributes(myc)\$args]))
cpar <- myc(val);
constrainpar <- c(constrainpar,list(cpar)); names(constrainpar) <- cname
if (p%in%ii\$parname.all) {
reg.idx <- which(x\$par==p)
reg.tidx <- which(t(x\$par==p))
if (!is.null(val))
A[reg.idx] <- cpar##myc(val)
}
if (p%in%ii\$covparname.all) {
cov.idx <- which(x\$covpar==p)
if (!is.null(val))
P[cov.idx] <- cpar##myc(val)
}
if (p%in%ii\$mparname.all) {
m.idx <- which(x\$mean==p)
if (!is.null(val))
v[m.idx] <- cpar##myc(val)
}
if (p%in%ii\$eparname.all) {
e.idx <- which(x\$exfix==p)
if (!is.null(val))
e[e.idx] <- cpar##myc(val)
}
constrain.idx[[p]] <- list(reg.idx=reg.idx,reg.tidx=reg.tidx,cov.idx=cov.idx,m.idx=m.idx,e.idx=e.idx)
}
}

if (x\$index\$sparse & !is.character(class(pars)[1])) {
A <- as(A,"sparseMatrix")
P <- as(P,"sparseMatrix")
v <- as(v,"sparseMatrix")
}
return(list(A=A, P=P, v=v, e=e, parval=parval, constrain.idx=constrain.idx, constrainpar=constrainpar))
}

###}}} matrices Obsolete
```
kkholst/lava documentation built on Sept. 6, 2021, 11:36 p.m.