mlnormal function

(Restricted) Maximum Likelihood Estimation with Prior Distributions and Penalty Functions under Multivariate Normality

The mlnormal estimates statistical model for multivariate normally distributed outcomes with specified mean structure and covariance structure (see Details and Examples). Model classes include multilevel models, factor analysis, structural equation models, multilevel structural equation models, social relations model and perhaps more.

The estimation can be conducted under maximum likelihood, restricted maximum likelihood and maximum posterior estimation with prior distribution. Regularization (i.e. LASSO penalties) is also accomodated.

mlnormal(y, X, id, Z_list, Z_index, beta=NULL, theta, method="ML", prior=NULL,
    lambda_beta=NULL, weights_beta=NULL, lambda_theta=NULL, weights_theta=NULL,
    beta_lower=NULL, beta_upper=NULL,    theta_lower=NULL, theta_upper=NULL,
    maxit=800, globconv=1e-05, conv=1e-06, verbose=TRUE, REML_shortcut=NULL,
    use_ginverse=FALSE, vcov=TRUE, variance_shortcut=TRUE, use_Rcpp=TRUE,
    level=0.95, numdiff.parm=1e-04, control_beta=NULL, control_theta=NULL)

## S3 method for class 'mlnormal'
summary(object, digits=4, file=NULL, ...)

## S3 method for class 'mlnormal'
print(x, digits=4, ...)

## S3 method for class 'mlnormal'
coef(object, ...)

## S3 method for class 'mlnormal'
logLik(object, ...)

## S3 method for class 'mlnormal'
vcov(object, ...)

## S3 method for class 'mlnormal'
confint(object, parm, level=.95, ... )

Arguments

• y: Vector of outcomes

• X: Matrix of covariates

• id: Vector of identifiers (subjects or clusters, see Details)

• Z_list: List of design matrices for covariance matrix (see Details)

• Z_index: Array containing loadings of design matrices (see Details). The dimensions are units $\times$ matrices $\times$ parameters.

• beta: Initial vector for $\bold{\beta}$

• theta: Initial vector for $\bold{\theta}$

• method: Estimation method. Can be either "ML" or "REML".

• prior: Prior distributions. Can be conveniently specified in a string which is processed by the function prior_model_parse. Only univariate prior distributions can be specified.

• lambda_beta: Parameter $\lambda_{\bold{\beta}}$ for penalty function $P( \bold{\beta} )=\lambda_{\bold{\beta}} \sum_h w_{\bold{\beta}h} | \beta _h |$

• weights_beta: Parameter vector $\bold{w}_{\bold{\beta}}$ for penalty function $P( \bold{\beta} )=\lambda_{\bold{\beta}} \sum_h w_{\bold{\beta}h} | \beta _h |$

• lambda_theta: Parameter $\lambda_{\bold{\theta}}$ for penalty function c("$P( \\bold{\\theta} )=\\lambda_{\\bold{\\theta}}\n$", "$\\sum_h w_{\\bold{\\theta}h} | \\theta _h |$")

• weights_theta: Parameter vector $\bold{w}_{\bold{\theta}}$ for penalty function c("$P( \\bold{\\theta} )=\\lambda_{\\bold{\\theta}}\n$", "$\\sum_h w_{\\bold{\\theta}h} | \\theta _h |$")

• beta_lower: Vector containing lower bounds for $\bold{\beta}$ parameter

• beta_upper: Vector containing upper bounds for $\bold{\beta}$ parameter

• theta_lower: Vector containing lower bounds for $\bold{\theta}$ parameter

• theta_upper: Vector containing upper bounds for $\bold{\theta}$ parameter

• maxit: Maximum number of iterations

• globconv: Convergence criterion deviance

• conv: Maximum parameter change

• verbose: Print progress?

• REML_shortcut: Logical indicating whether computational shortcuts should be used for REML estimation

• use_ginverse: Logical indicating whether a generalized inverse should be used

• vcov: Logical indicating whether a covariance matrix of $\bold{\theta}$ parameter estimates should be computed in case of REML (which is computationally demanding)

• variance_shortcut: Logical indicating whether computational shortcuts for calculating covariance matrices should be used

• use_Rcpp: Logical indicating whether the Rcpp package should be used

• level: Confidence level

• numdiff.parm: Numerical differentiation parameter

• control_beta: List with control arguments for $\bold{\beta}$ estimation. The default is

  list( maxiter=10, conv=1E-4, ridge=1E-6).

• control_theta: List with control arguments for $\bold{\theta}$ estimation. The default is

  list( maxiter=10, conv=1E-4, ridge=1E-6).

• object: Object of class mlnormal

• digits: Number of digits used for rounding

• file: File name

• parm: Parameter to be selected for confint method

• ...: Further arguments to be passed

• x: Object of class mlnormal

Details

The data consists of outcomes $\bold{y}_i$ and covariates $\bold{X}_i$

for unit $i$. The unit can be subjects, clusters (like schools) or the full outcome vector. It is assumed that $\bold{y}_i$ is normally distributed as $N( \bold{\mu}_i, \bold{V}_i )$ where the mean structure is modelled as

$$\bold{\mu}_i=\bold{X}_i \bold{\beta}$$

and the covariance structure $\bold{V}_i$ depends on a parameter vector $\bold{\theta}$. More specifically, the covariance matrix $\bold{V}_i$ is modelled as a sum of functions of the parameter $\bold{\theta}$ and known design matrices $\bold{Z}_{im}$ for unit $i$ ($m=1,\ldots,M$). The model is

$$\bold{V}_i=\sum_{m=1}^M \bold{Z}_{im} \gamma_{im} \qquad \mathrm{with}\qquad \gamma_{im}=\prod_{h=1}^H \theta_h^{q_{imh}}$$

where $q_{imh}$ are non-negative known integers specified in Z_index and $\bold{Z}_{im}$ are design matrices specified in Z_list.

The estimation follows Fisher scoring (Jiang, 2007; for applications see also Longford, 1987; Lee, 1990; Gill & Swartz, 2001) and the regularization approach is as described in Lin, Pang and Jiang (2013) (see also Krishnapuram, Carin, Figueiredo, & Hartemink, 2005).

Returns

List with entries

• theta: Estimated $\bold{\theta}$ parameter

• beta: Estimated $\bold{\beta}$ parameter

• theta_summary: Summary of $\bold{\theta}$ parameters

• beta_summary: Summary of $\bold{\beta}$ parameters

• coef: Estimated parameters

• vcov: Covariance matrix of estimated parameters

• ic: Information criteria

• V_list: List with fitted covariance matrices $\bold{V}_i$

• V1_list: List with inverses of fitted covariance matrices $\bold{V}_i$

• prior_args: Some arguments in case of prior distributions

• ...: More values

References

Gill, P. S., & Swartz, T. B. (2001). Statistical analyses for round robin interaction data. Canadian Journal of Statistics, 29, 321-331. tools:::Rd_expr_doi("10.2307/3316080")

Jiang, J. (2007). Linear and generalized linear mixed models and their applications. New York: Springer.

Krishnapuram, B., Carin, L., Figueiredo, M. A., & Hartemink, A. J. (2005). Sparse multinomial logistic regression: Fast algorithms and generalization bounds. IEEE Transactions on Pattern Analysis and Machine Intelligence, 27, 957-968. tools:::Rd_expr_doi("10.1109/TPAMI.2005.127")

Lee, S. Y. (1990). Multilevel analysis of structural equation models. Biometrika, 77, 763-772. tools:::Rd_expr_doi("10.1093/biomet/77.4.763")

Lin, B., Pang, Z., & Jiang, J. (2013). Fixed and random effects selection by REML and pathwise coordinate optimization. Journal of Computational and Graphical Statistics, 22, 341-355. tools:::Rd_expr_doi("10.1080/10618600.2012.681219")

Longford, N. T. (1987). A fast scoring algorithm for maximum likelihood estimation in unbalanced mixed models with nested random effects. Biometrika, 74, 817-827. tools:::Rd_expr_doi("10.1093/biomet/74.4.817")

See Also

See lavaan, sem, lava, OpenMx or nlsem

packages for estimation of (single level) structural equation models.

See the regsem

and lsl packages for regularized structural equation models.

See lme4 or nlme package for estimation of multilevel models.

See the lmmlasso and glmmLasso packages for regularized mixed effects models.

See OpenMx and xxM packages (http://xxm.times.uh.edu/) for estimation of multilevel structural equation models.

Examples

## Not run:

#############################################################################
# EXAMPLE 1: Two-level random intercept model
#############################################################################

#--------------------------------------------------------------
# Simulate data
#--------------------------------------------------------------

set.seed(976)
G <- 150 ; rg <- c(10,20)   # 150 groups with group sizes ranging from 10 to 20
#* simulate group sizes
ng <- round( stats::runif( G, min=rg[1], max=rg[2] ) )
idcluster <- rep(1:G, ng )
#* simulate covariate
iccx <- .3
x <- rep( stats::rnorm( G, sd=sqrt( iccx) ), ng ) +
            stats::rnorm( sum(ng), sd=sqrt( 1 - iccx) )
#* simulate outcome
b0 <- 1.5 ; b1 <- .4 ; iccy <- .2
y <- b0 + b1*x + rep( stats::rnorm( G, sd=sqrt( iccy) ), ng ) +
         stats::rnorm( sum(ng), sd=sqrt( 1 - iccy) )

#-----------------------
#--- arrange input for mlnormal function

id <- idcluster          # cluster is identifier
X <- cbind( 1, x )      # matrix of covariates
N <- length(id)          # number of units (clusters), which is G

MD <- max(ng)   # maximum number of persons in a group
NP <- 2         # number of covariance parameters theta

#* list of design matrix for covariance matrix
#  In the case of the random intercept model, the covariance structure is
#  tau^2 * J + sigma^2 * I, where J is a matrix of ones and I is the
#  identity matrix
Z <- as.list(1:G)
for (gg in 1:G){
    Ngg <- ng[gg]
    Z[[gg]] <- as.list( 1:2 )
    Z[[gg]][[1]] <- matrix( 1, nrow=Ngg, ncol=Ngg )  # level 2 variance
    Z[[gg]][[2]] <- diag(1,Ngg)            # level 1 variance
}
Z_list <- Z
#* parameter list containing the powers of parameters
Z_index <- array( 0, dim=c(G,2,2) )
Z_index[ 1:G, 1, 1] <- Z_index[ 1:G, 2, 2] <- 1

#** starting values and parameter names
beta <- c( 1, 0 )
names(beta) <- c("int", "x")
theta <- c( .05, 1 )
names(theta) <- c("tau2", "sig2" )

#** create dataset for lme4 for comparison
dat <- data.frame(y=y, x=x, id=id )

#--------------------------------------------------------------
# Model 1: Maximum likelihood estimation
#--------------------------------------------------------------

#** mlnormal function
mod1a <- LAM::mlnormal( y=y, X=X, id=id, Z_list=Z_list, Z_index=Z_index,
            beta=beta, theta=theta, method="ML" )
summary(mod1a)

# lme4::lmer function
library(lme4)
mod1b <- lme4::lmer( y ~ x + (1 | id ), data=dat, REML=FALSE )
summary(mod1b)

#--------------------------------------------------------------
# Model 2: Restricted maximum likelihood estimation
#--------------------------------------------------------------

#** mlnormal function
mod2a <- LAM::mlnormal( y=y, X=X, id=id, Z_list=Z_list, Z_index=Z_index,
            beta=beta, theta=theta, method="REML" )
summary(mod2a)

# lme4::lmer function
mod2b <- lme4::lmer( y ~ x + (1 | id ), data=dat, REML=TRUE )
summary(mod2b)

#--------------------------------------------------------------
# Model 3: Estimation of standard deviation instead of variances
#--------------------------------------------------------------

# The model is now parametrized in standard deviations
# Variances are then modeled as tau^2 and sigma^2, respectively.
Z_index2 <- 2*Z_index       # change loading matrix

# estimate model
mod3 <- LAM::mlnormal( y=y, X=X, id=id, Z_list=Z_list, Z_index=Z_index2,
            beta=beta, theta=theta )
summary(mod3)

#--------------------------------------------------------------
# Model 4: Maximum posterior estimation
#--------------------------------------------------------------

# specify prior distributions for parameters
prior <- "
    tau2 ~ dgamma(NA, 2, .5 )
    sig2 ~ dinvgamma(NA, .1, .1 )
    x ~ dnorm( NA, .2, 1000 )
    "

# estimate model in mlnormal
mod4 <- LAM::mlnormal( y=y, X=X, id=id, Z_list=Z_list, Z_index=Z_index,
            beta=beta, theta=theta, method="REML", prior=prior, vcov=FALSE )
summary(mod4)

#--------------------------------------------------------------
# Model 5: Estimation with regularization on beta and theta parameters
#--------------------------------------------------------------

#*** penalty on theta parameter
lambda_theta <- 10
weights_theta <- 1 + 0 * theta
#*** penalty on beta parameter
lambda_beta <- 3
weights_beta <- c( 0, 1.8 )

# estimate model
mod5 <- LAM::mlnormal( y=y, X=X, id=id, Z_list=Z_list, Z_index=Z_index,
            beta=beta, theta=theta, method="ML", maxit=maxit,
            lambda_theta=lambda_theta, weights_theta=weights_theta,
            lambda_beta=lambda_beta, weights_beta=weights_beta  )
summary(mod5)

#############################################################################
# EXAMPLE 2: Latent covariate model, two-level regression
#############################################################################

# Yb=beta_0 + beta_b*Xb + eb (between level) and
# Yw=beta_w*Xw + ew (within level)

#--------------------------------------------------------------
# Simulate data from latent covariate model
#--------------------------------------------------------------

set.seed(865)
# regression parameters
beta_0 <- 1 ; beta_b <- .7 ; beta_w <- .3
G <- 200      # number of groups
n <- 15      # group size
iccx <- .2   # intra class correlation x
iccy <- .35  # (conditional) intra class correlation y
# simulate latent variables
xb <- stats::rnorm(G, sd=sqrt( iccx ) )
yb <- beta_0 + beta_b * xb + stats::rnorm(G, sd=sqrt( iccy ) )
xw <- stats::rnorm(G*n, sd=sqrt( 1-iccx ) )
yw <- beta_w * xw + stats::rnorm(G*n, sd=sqrt( 1-iccy ) )
group <- rep( 1:G, each=n )
x <- xw + xb[ group ]
y <- yw + yb[ group ]
# test results on true data
lm( yb ~ xb )
lm( yw ~ xw )

# create vector of outcomes in the form
# ( y_11, x_11, y_21, x_21, ... )
dat <- cbind( y, x )
dat
Y <- as.vector( t(dat) )    # outcome vector
ny <- length(Y)
X <- matrix( 0, nrow=ny, ncol=2 )
X[ seq(1,ny,2), 1 ] <- 1   # design vector for mean y
X[ seq(2,ny,2), 2 ] <- 1   # design vector for mean x
id <- rep( group, each=2 )

#--------------------------------------------------------------
# Model 1: Linear regression ignoring multilevel structure
#--------------------------------------------------------------

# y=beta_0 + beta_t *x + e
# Var(y)=beta_t^2 * var_x + var_e
# Cov(y,x)=beta_t * var_x
# Var(x)=var_x

#** initial parameter values
theta <- c( 0, 1, .5 )
names(theta) <- c( "beta_t", "var_x", "var_e")
beta <- c(0,0)
names(beta) <- c("mu_y","mu_x")

# The unit i is a cluster in this example.

#--- define design matrices | list Z_list
Hlist <- list(  matrix( c(1,0,0,0), 2, 2 ), # var(y)
                matrix( c(1,0,0,0), 2, 2 ), # var(y) (two terms)
                matrix( c(0,1,1,0), 2, 2 ), # cov(x,y)
                matrix( c(0,0,0,1), 2, 2 ) ) # var(x)

U0 <- matrix( 0, nrow=2*n,ncol=2*n )
Ulist <- list( U0, U0, U0, U0 )
M <- length(Hlist)
for (mm in 1:M){    # mm <- 1
    for (nn in 1:n){     # nn <- 0
        Ulist[[ mm ]][ 2*(nn-1) + 1:2, 2*(nn-1) + 1:2 ] <- Hlist[[ mm ]]
    }
}
Z_list <- as.list(1:G)
for (gg in 1:G){
    Z_list[[gg]] <- Ulist
}

#--- define index vectors
Z_index <- array( 0, dim=c(G, 4, 3 ) )
K0 <- matrix( 0, nrow=4, ncol=3 )
colnames(K0) <- names(theta)
# Var(y)=beta_t^2 * var_x + var_e  (matrices withn indices 1 and 2)
K0[ 1, c("beta_t","var_x") ] <- c(2,1)  # beta_t^2 * var_x
K0[ 2, c("var_e") ] <- c(1)  # var_e
# Cov(y,x)=beta_t * var_x
K0[ 3, c("beta_t","var_x")] <- c(1,1)
# Var(x)=var_x
K0[ 4, c("var_x") ] <- c(1)
for (gg in 1:G){
    Z_index[gg,,] <- K0
}

#*** estimate model with mlnormal
mod1a <- LAM::mlnormal( y=Y, X=X, id=id, Z_list=Z_list, Z_index=Z_index,
            beta=beta, theta=theta, method="REML", vcov=FALSE )
summary(mod1a)

#*** estimate linear regression with stats::lm
mod1b <- stats::lm( y ~ x )
summary(mod1b)

#--------------------------------------------------------------
# Model 2: Latent covariate model
#--------------------------------------------------------------

#** initial parameters
theta <- c( 0.12, .6, .5, 0, .2, .2 )
names(theta) <- c( "beta_w", "var_xw", "var_ew",
                "beta_b", "var_xb", "var_eb")

#--- define design matrices | list Z_list
Hlist <- list(  matrix( c(1,0,0,0), 2, 2 ), # var(y)
                matrix( c(1,0,0,0), 2, 2 ), # var(y) (two terms)
                matrix( c(0,1,1,0), 2, 2 ), # cov(x,y)
                matrix( c(0,0,0,1), 2, 2 ) ) # var(x)
U0 <- matrix( 0, nrow=2*n,ncol=2*n )
Ulist <- list( U0, U0, U0, U0,  # within structure
               U0, U0, U0, U0  )  # between structure
M <- length(Hlist)
#*** within structure
design_within <- diag(n)  # design matrix within structure
for (mm in 1:M){    # mm <- 1
    Ulist[[ mm ]] <- base::kronecker( design_within, Hlist[[mm]] )
}
#*** between structure
design_between <- matrix(1, nrow=n, ncol=n)
      # matrix of ones corresponding to group size
for (mm in 1:M){    # mm <- 1
    Ulist[[ mm + M ]] <-  base::kronecker( design_between, Hlist[[ mm ]] )
}
Z_list <- as.list(1:G)
for (gg in 1:G){
    Z_list[[gg]] <- Ulist
}

#--- define index vectors Z_index
Z_index <- array( 0, dim=c(G, 8, 6 ) )
K0 <- matrix( 0, nrow=8, ncol=6 )
colnames(K0) <- names(theta)
# Var(y)=beta^2 * var_x + var_e  (matrices withn indices 1 and 2)
K0[ 1, c("beta_w","var_xw") ] <- c(2,1)  # beta_t^2 * var_x
K0[ 2, c("var_ew") ] <- c(1)  # var_e
K0[ 5, c("beta_b","var_xb") ] <- c(2,1)  # beta_t^2 * var_x
K0[ 6, c("var_eb") ] <- c(1)  # var_e
# Cov(y,x)=beta * var_x
K0[ 3, c("beta_w","var_xw")] <- c(1,1)
K0[ 7, c("beta_b","var_xb")] <- c(1,1)
# Var(x)=var_x
K0[ 4, c("var_xw") ] <- c(1)
K0[ 8, c("var_xb") ] <- c(1)
for (gg in 1:G){
    Z_index[gg,,] <- K0
}

#--- estimate model with mlnormal
mod2a <- LAM::mlnormal( y=Y, X=X, id=id, Z_list=Z_list, Z_index=Z_index,
            beta=beta, theta=theta, method="ML" )
summary(mod2a)

#############################################################################
# EXAMPLE 3: Simple linear regression, single level
#############################################################################

#--------------------------------------------------------------
# Simulate data
#--------------------------------------------------------------

set.seed(875)
N <- 300
x <- stats::rnorm( N, sd=1.3 )
y <- .4 + .7 * x + stats::rnorm( N, sd=.5 )
dat <- data.frame( x, y )

#--------------------------------------------------------------
# Model 1: Linear regression modelled with residual covariance structure
#--------------------------------------------------------------

# matrix of predictros
X <- cbind( 1, x )
# list with design matrices
Z <- as.list(1:N)
for (nn in 1:N){
    Z[[nn]] <- as.list( 1 )
    Z[[nn]][[1]] <- matrix( 1, nrow=1, ncol=1 )  # residual variance
}
#* loading matrix
Z_index <- array( 0, dim=c(N,1,1) )
Z_index[ 1:N, 1, 1] <- 2  # parametrize residual standard deviation
#** starting values and parameter names
beta <- c( 0, 0 )
names(beta) <- c("int", "x")
theta <- c(1)
names(theta) <- c("sig2" )
# id vector
id <- 1:N

#** mlnormal function
mod1a <- LAM::mlnormal( y=y, X=X, id=id, Z_list=Z, Z_index=Z_index,
            beta=beta, theta=theta, method="ML" )
summary(mod1a)

# estimate linear regression with stats::lm
mod1b <- stats::lm( y ~ x )
summary(mod1b)

#--------------------------------------------------------------
# Model 2: Linear regression modelled with bivariate covariance structure
#--------------------------------------------------------------

#** define design matrix referring to mean structure
X <- matrix( 0, nrow=2*N, ncol=2 )
X[ seq(1,2*N,2), 1 ] <- X[ seq(2,2*N,2), 2 ] <- 1

#** create outcome vector
y1 <- dat[ cbind( rep(1:N, each=2), rep(1:2, N ) ) ]
#** list with design matrices
Z <- as.list(1:N)
Z0 <- 0*matrix( 0, nrow=2,ncol=2)
ZXY <- ZY <- ZX <- Z0
# design matrix Var(X)
ZX[1,1] <- 1
# design matrix Var(Y)
ZY[2,2] <- 1
# design matrix covariance
ZXY[1,2] <- ZXY[2,1] <- 1
# Var(X)=sigx^2
# Cov(X,Y)=beta * sigx^2
# Var(Y)=beta^2 * sigx^2 + sige^2
Z_list0 <- list( ZY, ZY, ZXY, ZX )
for (nn in 1:N){
    Z[[nn]] <- Z_list0
}
#* parameter list containing the powers of parameters
theta <- c(1,0.3,1)
names(theta) <- c("sigx", "beta", "sige" )
Z_index <- array( 0, dim=c(N,4,3) )
for (nn in 1:N){
    # Var(X)
    Z_index[nn, 4, ] <- c(2,0,0)
    # Cov(X,Y)
    Z_index[nn, 3, ] <- c(2,1,0)
    # Var(Y)
    Z_index[nn,1,] <- c(2,2,0)
    Z_index[nn,2,] <- c(0,0,2)
}
#** starting values and parameter names
beta <- c( 0, 0 )
names(beta) <- c("Mx", "My")
# id vector
id <- rep( 1:N, each=2 )

#** mlnormal function
mod2a <- LAM::mlnormal( y=y1, X=X, id=id, Z_list=Z, Z_index=Z_index,
            beta=beta, theta=theta, method="ML" )
summary(mod2a)

#--------------------------------------------------------------
# Model 3: Bivariate normal distribution in (sigma_X, sigma_Y, sigma_XY) parameters
#--------------------------------------------------------------

# list with design matrices
Z <- as.list(1:N)
Z0 <- 0*matrix( 0, nrow=2,ncol=2)
ZXY <- ZY <- ZX <- Z0
# design matrix Var(X)
ZX[1,1] <- 1
# design matrix Var(Y)
ZY[2,2] <- 1
# design matrix covariance
ZXY[1,2] <- ZXY[2,1] <- 1
Z_list0 <- list( ZX, ZY, ZXY  )
for (nn in 1:N){
    Z[[nn]] <- Z_list0
}

#* parameter list
theta <- c(1,1,.3)
names(theta) <- c("sigx", "sigy", "sigxy" )
Z_index <- array( 0, dim=c(N,3,3) )
for (nn in 1:N){
    # Var(X)
    Z_index[nn, 1, ] <- c(2,0,0)
    # Var(Y)
    Z_index[nn, 2, ] <- c(0,2,0)
    # Cov(X,Y)
    Z_index[nn, 3, ] <- c(0,0,1)
}

#** starting values and parameter names
beta <- c( 0, 0 )
names(beta) <- c("Mx", "My")

#** mlnormal function
mod3a <- LAM::mlnormal( y=y1, X=X, id=id, Z_list=Z, Z_index=Z_index,
            beta=beta, theta=theta, method="ML" )
summary(mod3a)

#--------------------------------------------------------------
# Model 4: Bivariate normal distribution in parameters of Cholesky decomposition
#--------------------------------------------------------------

# list with design matrices
Z <- as.list(1:N)
Z0 <- 0*matrix( 0, nrow=2,ncol=2)
ZXY <- ZY <- ZX <- Z0
# design matrix Var(X)
ZX[1,1] <- 1
# design matrix Var(Y)
ZY[2,2] <- 1
# design matrix covariance
ZXY[1,2] <- ZXY[2,1] <- 1
Z_list0 <- list( ZX, ZXY, ZY, ZY  )
for (nn in 1:N){
    Z[[nn]] <- Z_list0
}

#* parameter list containing the powers of parameters
theta <- c(1,0.3,1)
names(theta) <- c("L11", "L21", "L22" )
Z_index <- array( 0, dim=c(N,4,3) )
for (nn in 1:N){
    Z_index[nn,1,] <- c(2,0,0)
    Z_index[nn,2,] <- c(1,1,0)
    Z_index[nn,3,] <- c(0,2,0)
    Z_index[nn,4,] <- c(0,0,2)
}
#** starting values and parameter names
beta <- c( 0, 0 )
names(beta) <- c("Mx", "My")
# id vector
id <- rep( 1:N, each=2 )
#** mlnormal function
mod4a <- LAM::mlnormal( y=y1, X=X, id=id, Z_list=Z, Z_index=Z_index,
            beta=beta, theta=theta, method="ML" )
# parameter with lower diagonal entries of Cholesky matrix
mod4a$theta
# fill-in parameters for Cholesky matrix
L <- matrix(0,2,2)
L[ ! upper.tri(L) ] <- mod4a$theta
#** reconstruct covariance matrix
L 
stats::cov.wt(dat, method="ML")$cov
## End(Not run)