R/within_varcov.r

In mars: Meta Analysis and Research Synthesis

within_varcov <- function(data, N, effect_name, study, type,
                          variance = NULL, ...) {

  if(!type %in% c('outcome', 'time', 'control', 'multilevel',
                  'univariate')) {
    data_list <- split(data, f = data[[study]])
    cor_matrix <- lapply(seq_along(data_list), function(xx)
      corpcor::vec2sm(data_list[[xx]][[effect_name]]))
    for(xx in seq_along(cor_matrix)) {
      diag(cor_matrix[[xx]]) <- 1
    }
  }

  switch(type,
        simple = olkin_simple(cor_matrix, N),
        weighted = olkin_weighted(cor_matrix, N),
        average = olkin_average(cor_matrix, N),
        outcome = smd_outcome(data, effect_name, study),
        time = smd_time(data),
        control = smd_control(data),
        multilevel = multilevel_var(data,
                                    variance,
                                    study),
        univariate = univariate_var(data,
                                    variance,
                                    study)
        )
}

univariate_var <- function(data, variance, study) {
  S_list <- split(data[c(variance, study)], f = data[study])
  lapply(seq_along(S_list), function(xx)
    S_list[[xx]][[variance]])
}

multilevel_var <- function(data, variance,
                           study) {
  S_list <- split(data[c(variance, study)], f = data[study])

  lapply(seq_along(S_list), function(xx)
    diag(S_list[[xx]][[variance]]))

}

olkin_simple <- function(data, N) {
  olkin_siotani(data, N, type = 'simple')
}

olkin_weighted <- function(data, N) {
  olkin_siotani(data, N, type = 'weighted')
}

olkin_average <- function(data, N) {
  olkin_siotani(data, N, type = 'average')
}

smd_time <- function(data) {

}

smd_outcome <- function(data, effect_name, study) {
  WithinVarCov(data, effect_name, study)
}

smd_control <- function(data) {

}

datPairs <-  function(data, effect_name) {
  allPair <- combn(data[[effect_name]],2)
  allPairData <- data.frame(cbind(t(allPair),
                                  rep(data[['nt']][1], length(allPair[1,])),
                                  rep(data[['nc']][1], length(allPair[1,])),
                                  rep(data[['r2do']][1], length(allPair[1,]))))
  names(allPairData) <- c('dt', 'dc', 'nt', 'nc', 'r2do')
  return(allPairData)
}

#----- multiple ds outco ----------------------
varD <- function(data, effect_name) { # variance within
  myvar <- ( (1/data['nc']) + (1/data[,'nt'])) +
    data[[effect_name]]^2 / (2*( (data['nc']) + (data['nt'])))
  return(as.numeric(myvar))
}

dCov2Out <- function(data){ # covariance within
  mycov <- (( (1/data['nc']) + (1/data['nt'])) * data['r2do'] ) +
    ( (.5 *data['dt'] * data['dc'] * data['r2do']^2 ) /
        sum(data['nc'], data['nt']) )
  return(as.numeric(mycov))
}

WithinVarCov <- function(data, effect_name, study) {

  myList <- split(data, data[[study]])

  #----------- building mat -----------------
  withinVar <-  lapply(myList, function(x){
    res <-  NULL
    for (i in 1:nrow(x)) {
      res = c(res, varD(x[i,], effect_name))
    }
    return(res)}
  )

  withinCov <- lapply(lapply(myList, datPairs, effect_name = effect_name), function(x) apply(x, 1, dCov2Out))
  #----------------- computing var cov -----------------------
  res2 <- vector("list", length(myList))
  for (j in 1:length(myList)) {
    dat1 <-  withinCov[[j]]
    dat2 <- withinVar[[j]]
    res2[[j]] <- corpcor::vec2sm(dat1, diag = FALSE, order = NULL)
    diag(res2[[j]]) <- dat2
  }

  return(res2)
}