Suppose we look at the following model
$$ \vec y_i=\vec\mu_i + \vec\epsilon_i, \qquad \vec\epsilon_i\sim N(\vec0, \Sigma) $$
where $\vec y_i$s is observed, $\vec\mu_i$s are known, and $\vec\epsilon_i$s are unknown and iid. Following
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-PLUS. New York: Springer.
then we can we can separate $\Sigma$ into standard deviation factor and correlation matrix factor
$$ \Sigma =VCV, \quad V = \begin{pmatrix} \sigma_1 & 0 & \cdots & 0 \\ 0 & \sigma_2 & \ddots & 0 \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & \sigma_p \end{pmatrix}, \quad C = \begin{pmatrix} 1 & \rho_{21} & \cdots & \rho_{p1} \\ \rho_{21} & 1 & \ddots & \rho_{p2} \\ \vdots & \ddots & \ddots & \vdots \\ \rho_{p1} & \dots & \rho_{p,p-1} & 1 \end{pmatrix} $$
and restrict the $p(p+1)/2$ parameters above. In particular, say that we model
$$ \begin{align*} \sigma_{i} &= \exp(s_i), \quad (s_1,s_2,\dots, s_p)^\top=K\vec\psi \\ \rho_{ij}&= \frac{2}{1 + \exp(-q_{ij})} - 1, \quad (q_{21},q_{31},\dots,q_{p1},q_{32},\dots,q_{p,p-1})^\top=L\vec\varphi \end{align*}$$
where $K\in\mathbb{R}^{p\times k}$, $L\in\mathbb{R}^{p(p-1)/2\times l}$, both are known, $k\leq p$ and $l\leq p(p-1)/2$. An issue with this approach is that $C$ may not be correlation matrix. However, computing the log-likelihood is straightforward and so is computing the gradient. Also it allows for quite general covariance matrices.
On the flip site, one can implement a series of different model as in the nlme package some of which are not included in the above. This will take time but one can write better optimized function for each specific case.
My questions are:
- Is the above a bad idea?
- Is the above expected to often failed when using numerical optimizer? Can one do better in a general model?
- Is there a smarter and more general way to restrict the parameters and still have quite general models?
I am considering the above to allow the end user to restrict parameters in an R package (which is not pure multivariate normal model as above) and the above seems easy for the end user to specify. An obvious assumption of the end user would be to restrict some of the correlations to be zero. References will be appreciated.
Update
It is not clear from the above in what context I use the model. It is used in a hidden Markov model as part of a Monte Carlo expectation conditional maximization algorithm. See the example here where I have implemented the above in one of the conditional maximization steps. $l$ and $k$ are typically small and $p$ may be small or moderate.
Here is an example to make it more concrete. First I show the output and then the two function definitions (get_covar and ex_func)
######
# first example
K <- matrix(1, 4, 1)
psi <- log(2)
L <- matrix(0, 4 * (4 - 1) / 2, 1)
varphi <- log(- (.5 + 1) / (.5 - 1))
L[2, 1] <- 1
get_covar(K, L, psi, varphi)$Q
#R [,1] [,2] [,3] [,4]
#R [1,] 4 0 2 0
#R [2,] 0 4 0 0
#R [3,] 2 0 4 0
#R [4,] 0 0 0 4
set.seed(69832532)
out <- ex_func(K, L, psi, varphi)
do.call(rbind, lapply(out, "[[", "par"))
#R [,1] [,2]
#R wo_gr 0.6849 1.025
#R w_gr 0.6849 1.025
c(psi, varphi) # actual values
#R [1] 0.6931 1.0986
#####
# second example
K <- matrix(0, 4, 2)
K[1:2, 1] <- K[3:4, 2] <- 1
psi <- log(c(2, 5))
L <- matrix(0, 4 * (4 - 1) / 2, 2)
L[c(1, 2, 4), 1] <- 1
L[6, 2] <- 1
varphi <- log(- (c(.8, .4) + 1) / (c(.8, .4) - 1))
get_covar(K, L, psi, varphi)$Q
#R [,1] [,2] [,3] [,4]
#R [1,] 4.0 3.2 8 0
#R [2,] 3.2 4.0 8 0
#R [3,] 8.0 8.0 25 10
#R [4,] 0.0 0.0 10 25
set.seed(93900343)
out <- ex_func(K, L, psi, varphi)
do.call(rbind, lapply(out, "[[", "par"))
#R [,1] [,2] [,3] [,4]
#R wo_gr 0.7069 1.613 2.154 0.8533
#R w_gr 0.7070 1.613 2.154 0.8533
c(psi, varphi) # actual values
#R [1] 0.6931 1.6094 2.1972 0.8473
Here are the two function definitions
# function to get the covariance matrix
get_covar <- function(K, L, psi, varphi){
V <- diag(exp(drop(K %*% psi)))
C <- diag(1, ncol(V))
C[lower.tri(C)] <- 2/(1 + exp(-drop(L %*% varphi))) - 1
C[upper.tri(C)] <- t(C)[upper.tri(C)]
list(Q = V %*% C %*% V, V = V, C = C)
}
# function to simulate data and find MLE estimates
ex_func <- function(K, L, psi, varphi, nobs = 100){
# get covariance matrix
Q <- get_covar(K, L, psi, varphi)$Q
p <- ncol(Q)
# simulate some mus... Though, does not matter
mus <- matrix(rnorm(nobs * p), ncol = p)
Y <- mus + matrix(rnorm(nobs * p), ncol = p) %*% chol(Q)
# ... since we do
Z <- crossprod(Y - mus)
# assign log-likelihood function
ll <- function(par, k = length(psi)){
idx <- 1:k
psi <- par[ idx]
varphi <- par[-idx]
Q <- get_covar(K, L, psi, varphi)$Q
Q_qr <- qr(Q)
deter <- determinant(Q, logarithm = TRUE)
if(deter$sign < 0 || Q_qr$rank < ncol(Q))
return(NA_real_)
-(nobs * deter$modulus + sum(diag(solve(Q_qr, Z)))) / 2
}
# assign gradient function
gr <- function(par, k = length(psi)){
idx <- 1:k
psi <- par[ idx]
varphi <- par[-idx]
tmp <- get_covar(K, L, psi, varphi)
list2env(tmp, environment())
# Computations could be done a lot smarter...
Q_qr <- qr(Q)
if(Q_qr$rank < ncol(Q))
return(NA_real_)
fac <- solve(Q_qr, Z) - diag(nobs, ncol(Q))
fac <- solve(Q_qr, t(fac)) / 2
d_V <-
diag(fac %*% V %*% C + C %*% V %*% fac) %*%
K %*% diag(exp(psi), length(psi))
d_C <- tcrossprod(V, V %*% fac)
exp_varphi <- exp(varphi)
d_C <- as.vector(d_C)[lower.tri(d_C, diag = FALSE)] %*% L %*%
diag((4 * exp_varphi / (1 + exp_varphi)^2), length(varphi))
c(drop(d_V), drop(d_C))
}
# use optim
# scale with full covariance matrix log-likelihood
Q_full <- Z / nobs
deter <- determinant(Q_full, logarithm = TRUE)
ll_full <- -(nobs * deter$modulus + nobs * ncol(Q_full)) / 2
ctrl <- list(fnscale = if(ll_full < 0) ll_full else -ll_full)
start <- c(rep(1, length(psi)), rep(0, length(varphi)))
list(
wo_gr = optim(start, ll, control = ctrl ),
w_gr = optim(start, ll, gr, control = ctrl, method = "BFGS"))
}
