R/model_stats.r

In mars: Meta Analysis and Research Synthesis

b_fe <- function(X, invS, Y) {
  chol2inv(chol(t(X) %*% invS %*% X)) %*% t(X) %*% invS %*% Y
}


QE <- function(X, b, Y, invS) {
  t(Y - X %*% b) %*% invS %*% (Y - X %*% b)
}


QE_compute <- function(X, b, Y, invS) {
  QE_value <- QE(X, b, Y, invS)
  QE_df <- length(Y) - length(b)
  QE_p <- pchisq(q = QE_value, df = QE_df, lower.tail = FALSE)

  # Brainstorm this
  list(value = QE_value,
       df = QE_df,
       p = QE_p)
}


QE_within <- function(data, effectID, effectsize_name, S_name) {
  data_list <- split(data, data[effectID])
  QE_out <- lapply(data_list, function(xx) {
    if( nrow(xx) > 1){
      si <- xx[, S_name]
      S <- diag(xx[, S_name])
      Y <- xx[[effectsize_name]]
      X <- matrix(rep(1, nrow(xx)))
      b <- chol2inv(chol(t(X) %*% chol2inv(chol(S)) %*% X)) %*% t(X) %*% chol2inv(chol(S)) %*% Y
      wi <- 1/si
      QE_value <- sum(wi * ((Y - rep(b, nrow(xx))))^2)
      QE_df <- length(Y) - 1
      QE_p <- pchisq(q = QE_value, df = QE_df, lower.tail = FALSE)
      c("QE_value" = QE_value, "QE_df" = QE_df, "QE_p" = QE_p)
    } else { "only one effect size in this category" }
  }
  )
  return(QE_out)
}



# I2

# multivariate for each outcome using common average
I2_within <- function(X, invS, Tau_est, q) {
  P <- invS - invS %*% X %*% chol2inv(chol(t(X) %*% invS %*% X)) %*% t(X) %*% invS
  100 * diag(Tau_est) / (diag(Tau_est) + (nrow(X) - q)/sum(diag(P)))
}


# Jackson et al. (2012)
I2_jackson <- function(X, invS, mu_cov) {
  100 * (diag(mu_cov) - diag(chol2inv(chol(t(X) %*% invS %*% X))) )/ diag(mu_cov)
}

# Overall I2 from average between
I2_between <- function(X, invS, Tau_est, q) {
  P <- invS - invS %*% X %*% chol2inv(chol(t(X) %*% invS %*% X)) %*% t(X) %*% invS
  100 * sum(diag(Tau_est))/q / ( sum(diag(Tau_est))/q + (nrow(X) - q)/sum(diag(P)))
}

# QM

QM <- function(mu_est, mu_cov) {
  QM_value <- t(matrix(mu_est)) %*% chol2inv(chol(mu_cov)) %*% mu_est
  QM_df <- length(mu_est)
  QM_p <- pchisq(q = QM_value, df = QM_df, lower.tail = FALSE)

  list(value = QM_value,
       df = QM_df,
       p = QM_p)
}



# need to get the weights with the mixed effects
# Tau.est
# Ss
inverse_psi <- function(varcov, Tau_est, num_studies, yID, Z_list = NULL,
                        structure, inverse = TRUE) {

  Tau <- diag(Tau_est)
  if(!is.null(Z_list)) {
      Psi <- lapply(seq_along(Z_list), function(i)
        sum_list(lapply(seq_along(Tau),function(j)
          as.matrix(Matrix::bdiag(lapply(Z_list[[i]][[j]], function(x)
            x %*% Tau[[j]] %*% t(x)))))) + varcov[[i]])
  } else {
    if(is.null(Tau_est)) {
      Psi <- vector(mode = "list", length = num_studies)
      if(structure == 'univariate') {
        for(i in seq_len(num_studies)){
          Psi[[i]] <- varcov[[i]][yID[[i]]]
        }
      } else {
        for(i in seq_len(num_studies)){
          Psi[[i]] <- varcov[[i]][yID[[i]], yID[[i]]]
        }
      }
    } else {
      if(structure == 'univariate') {
        Psi <- vector(mode = "list", length = num_studies)
        for(i in seq_len(num_studies)){
          Psi[[i]] <- varcov[[i]][yID[[i]]]  + Tau_est[yID[[i]]]
        }
      } else {
        Psi <- vector(mode = "list", length = num_studies)
        for(i in seq_len(num_studies)){
          Psi[[i]] <- varcov[[i]][yID[[i]], yID[[i]]]  + Tau_est[yID[[i]], yID[[i]] ]
        }
      }
    }
  }

  if(inverse) {
    if(structure == 'univariate') {
      chol2inv(chol(diag(unlist(Psi), nrow = num_studies)))
    } else {
      chol2inv(chol(Matrix::bdiag(Psi)))
    }
  } else {
    if(structure == 'univariate') {
      diag(unlist(Psi), nrow = num_studies)
    } else {
      Matrix::bdiag(Psi)
    }

  }

}



b_mix <- function(X, invPsi, Y) { # with intercept
  if(!is.matrix(X)) {
    X <- as.matrix(X)
  }
  if(!is.matrix(invPsi)) {
    invPsi <- as.matrix(invPsi)
  }

  part1 <- t(X) %*% invPsi %*% X

  # if(!matrixcalc::is.positive.definite(part1)) {
  #   part1 <- Matrix::nearPD(part1)[['mat']]
  # }
  chol2inv(chol(part1)) %*% t(X) %*% invPsi %*% Y
}


covb_mix <- function(X, invPsi) { # with intercept
  chol2inv(chol(t(X) %*% invPsi %*% X))
}

# need to compute h2 from
QE_mix <- function(X, bmix, Y, invPsi){ # only for the Dtest
  t(Y - X %*% bmix) %*% invPsi %*% (Y - X %*% bmix)
}

#' @importFrom stats pf
f_test <- function(QMmix, QEmix, num_obs) {
  # QMmix <- t(matrix(bmix)) %*% chol2inv(chol(covbmix)) %*% bmix
  # QEmix <- QE_mix(X, bmix, Y, invPsi)

  F_value <- QMmix[['value']] / QEmix

  df1 <- QMmix[['df']]
  df2 <- num_obs - QMmix[['df']]
  Fp_value <- pf(F_value, df1, df2, lower.tail = FALSE)

  list(value = F_value,
       df1 = df1,
       df2 = df2,
       p = Fp_value)
}

predicted_values <- function(beta, design_matrix) {
  design_matrix %*% beta
}

residual_values <- function(predicted_values, effect_size) {
  effect_size - predicted_values
}

robust_se <- function(data, clusterID, Psi, X,
                      cov_b, residuals) {

  clstr <- data[, clusterID]
  indctr <- order(clstr)
  clstr <- clstr[indctr]


  ### if no werror_ights were specified, then vb = (X'WX)^-1, so we can use that part

  #W <- chol2inv(chol(res$M[indctr,indctr]))   # inverse of the block diagonal
  invPsi_rve <- chol2inv(chol(Psi[indctr, indctr]))
  #bread <- res$vb %*% crossprod(res$X[indctr,], W) # var-coc slopes times cross product of indicator matrix and werror_ights

  bread <- cov_b %*% crossprod(X[indctr,], invPsi_rve)



  ### construct meat part

  error_i <- residuals[indctr]

  clstr <- factor(clstr, levels = unique(clstr))

  meat <- Matrix::bdiag(lapply(split(error_i, clstr), function(xx) tcrossprod(xx)))

  ###  robust var-cov matrix
  vb <- bread %*% meat %*% t(bread)

  n <- length(unique(clstr))

  dfs <- n - length(ncol(X)) # columns of design matrix

  vb_1 <- (n / ( n - 1)) * vb
  vb_2 <- ((n / (n - 1)) * ((nrow(X) - 1) / (nrow(X) - length(ncol(X))))  ) * vb

  list(
    vb_1 = vb_1,
    vb_2 = vb_2
  )

}


#---------- for CI --------------------------

#' @importFrom stats qnorm
ciFun <- function(beta, level, se) { # two for each one lower one upper
  crit <- qnorm(level/2, lower.tail=FALSE)
  Low <- c(beta - crit * se)
  Upp <- c(beta + crit * se)
  c(Low, Upp)
}

#ciFun(b, .95, sqrt(diag(covbmix)))