R/correlation_structures.r

In mars: Meta Analysis and Research Synthesis

correlation_structures <- function(params, q_r, q_f, structure = 'UN') {
  switch(structure,
         UN = unstructured(params, q_r = q_r),
         DIAG1 = diag1(params, q_r = q_r, q_f = q_f),
         DIAG2 = diag2(params, q_r = q_r, q_f = q_f),
         CS = compound_symmetry(params, q_r = q_r, q_f = q_f),
         HCS = h_compound_symmetry(params, q_r = q_r, q_f = q_f),
         AR1 = ar1(params, q_r = q_r, q_f = q_f),
         univariate = diag1(params, q_r = q_r, q_f = q_f)
  )
}

unstructured <- function(params, q_r) {

  L_tau <- matrix(0, nrow = q_r, ncol = q_r)
  for(i in 1:q_r){
    for(j in 1:i){
      L_tau[i, j] <- params[ q_r + (i - 1) * i/2 + j]
    }
  }
  tcrossprod(L_tau)

}

diag1 <- function(params, q_r, q_f) {

  diag(params[-1:-q_f], q_r)

}

diag2 <- function(params, q_r, q_f) {
  diag(rep(params[-1:-q_f], q_r), q_r)
}

compound_symmetry <- function(params, q_r, q_f) {
  R <- matrix(params[q_f + 2], q_r, q_r)
  diag(R) <- 1
  params[q_f + 1]  * R
}

h_compound_symmetry <- function(params, q_r, q_f) {
  R <- matrix(rep(params[q_f + q_r + 1]), q_r, q_r)
  diag(R) <- 1
  D <- diag(sqrt(params[q_f + seq(q_r)]*2), q_r)
  D %*% R %*% D
}

ar1 <- function(params, q_r, q_f) {
  expo <- abs(matrix(1:q_r - 1, q_r, q_r, byrow = TRUE) - (1:q_r - 1))
  R <- params[q_f + 2]^expo
  D <- diag(sqrt(params[q_f + 1]), q_r)
  D %*% R %*% D
}


correlation_export <- function(estimated_pars,
                               q_f = NULL,
                               q_r = NULL,
                               structure = 'UN') {
  switch(structure,
         UN = unstructured_est(estimated_pars, q_r = q_r),
         DIAG1 = diag1_est(estimated_pars, q_f = q_f, q_r = q_r),
         DIAG2 = diag2_est(estimated_pars, q_f = q_f, q_r = q_r),
         CS = compound_symmetry_est(estimated_pars, q_f = q_f, q_r = q_r),
         HCS = h_compound_symmetry_est(estimated_pars, q_f = q_f, q_r = q_r),
         AR1 = ar1_est(estimated_pars, q_f = q_f, q_r = q_r),
         multilevel = multilevel_est(estimated_pars, q_f = q_f, q_r = q_r),
         univariate = multilevel_est(estimated_pars, q_f = q_f, q_r = q_r)
  )
}

unstructured_est <- function(estimated_pars, q_r) {

  L_tau_est <- matrix(0, nrow = q_r, ncol = q_r)
  for(i in 1:q_r) {
    for(j in 1:i) {
      L_tau_est[i, j] <- estimated_pars[q_r + (i - 1)*i/2 + j]
    }
  }

  tcrossprod(L_tau_est)

}

diag1_est <- function(estimated_pars, q_f, q_r) {
  diag(estimated_pars[q_f + 1:q_r], q_r)
}

diag2_est <- function(estimated_pars, q_f, q_r) {
  diag(estimated_pars[q_f + 1], q_r)
}

multilevel_est <- function(estimated_pars, q_f, q_r) {
  diag(estimated_pars[(q_f+1):(q_f+q_r)], q_r)
}

compound_symmetry_est <- function(estimated_pars, q_f, q_r) {
  Tau <- matrix(estimated_pars[q_f + 2] * estimated_pars[q_f + 1] * estimated_pars[q_f + 1],
                    nrow = q_r, ncol = q_r)
  diag(Tau) <- estimated_pars[q_f + 1]
  Tau
}

h_compound_symmetry_est <- function(estimated_pars, q_f, q_r) {
  R <- matrix(rep(estimated_pars[q_f + q_r + 1]), q_r, q_r)
  diag(R) <- 1
  D <- diag(sqrt(estimated_pars[q_f + seq(q_r)]*2), q_r)
  D %*% R %*% D
}

ar1_est <- function(estimated_pars, q_f, q_r) {

  expo <- abs(matrix(1:q_r - 1, q_r, q_r, byrow = TRUE) - (1:q_r - 1))
  R <- estimated_pars[q_f + 2]^expo
  D <- diag(sqrt(estimated_pars[q_f + 1]), q_r)
  D %*% R %*% D
}