3
$\begingroup$

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"))
}
$\endgroup$
1
$\begingroup$

I think it is a very interesting idea but very difficult to implement, due to the positive definiteness constraint. What dimensionality are we talking about for $y$ and $\phi$?

I didn't recognise the formula for $\rho_{ij}$ as the inverse of the Fisher transform. Is this what it is intended to be?

What type of data are you modelling? If it is longitudinal, then there is another approach based on the modified Cholesky decomposition, implemented in the R package jmcm.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ $\vec{y}_i$ may be low or moderate dimensional and $l$ and $k$ are small. See the updated questions. No it was a typo with $q$. Thanks for correcting it. It is in a hidden Markov model. I will checkout the JSS paper for the package you mention. $\endgroup$ – Benjamin Christoffersen Aug 15 '18 at 9:50
  • $\begingroup$ I meant the formula for $\rho_{ij}$, not $q$. $\endgroup$ – papgeo Aug 15 '18 at 10:08
  • $\begingroup$ They are a function of $q_{ij}$ which are a function of $\vec{\varphi}$ as I write. $\endgroup$ – Benjamin Christoffersen Aug 15 '18 at 10:30
  • $\begingroup$ Section "2.4. Hyperspherical parameterization of the Cholesky factor (HPC)" in the JSS paper for the jmcmc package seems like one approach. However, it is not easy for me to see how one can imply restrictions on the correlation matrix $R_i$ elements through $\theta_{ijk}$. E.g., specifying that $R_{ijk} = 0$ seems complicated. Maybe I am wrong. $\endgroup$ – Benjamin Christoffersen Aug 15 '18 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.