library(RSiena)
library(metafor)
library(xlsx)


# rename result objects from selected model specification
for (i in village.codes) {
  
  objname <- paste0('saom.', i, '.rate', rates, '.', saom.model)
  assign(paste0('saom.result.', i), eval(parse(text=objname)))
  
}
rm(i, objname)

# set up table for meta-analysis results
saom.covars <- saom.result.v1$effects$effectName
saom.meta <- matrix(NA, length(saom.covars), 7)
rownames(saom.meta) <- saom.covars
colnames(saom.meta) <- c('mu', 'se', 'p', 'tau', 'Q', 'Qp', 'n')

# perform meta-analysis
for(i in 1:nrow(saom.meta)) {
  if (rates<8) {
    estimates <- c(saom.result.v1$theta[i], saom.result.v2$theta[i], saom.result.v3$theta[i],
                   saom.result.v4$theta[i], saom.result.v5$theta[i], saom.result.v6$theta[i],
                   saom.result.v7$theta[i], saom.result.v8$theta[i], saom.result.v9$theta[i],
                   saom.result.v10$theta[i])
    serrors <- c(saom.result.v1$se[i], saom.result.v2$se[i], saom.result.v3$se[i],
                 saom.result.v4$se[i], saom.result.v5$se[i], saom.result.v6$se[i],
                 saom.result.v7$se[i], saom.result.v8$se[i], saom.result.v9$se[i],
                 saom.result.v10$se[i])
  } else {
    estimates <- c(saom.result.v1$theta[i], saom.result.v2$theta[i], saom.result.v3$theta[i],
                   saom.result.v5$theta[i])
    serrors <- c(saom.result.v1$se[i], saom.result.v2$se[i], saom.result.v3$se[i],
                 saom.result.v5$se[i])
  }
  if (sum(is.na(serrors))<10) {
    meta <- rma(estimates, serrors^2)
    saom.meta[i,'mu'] <- meta$b
    saom.meta[i,'se'] <- meta$se
    saom.meta[i,'p'] <- meta$pval
    saom.meta[i,'tau'] <- sqrt(meta$tau2)
    saom.meta[i,'Q'] <- meta$QE
    saom.meta[i,'Qp'] <- meta$QEp
    saom.meta[i,'n'] <- meta$k
  }
}
rm(i, estimates, serrors, meta)

# rename meta-analysis object to reflect rate parameters and specification
objname <- paste0('saom.meta.rate', rates, '.', saom.model)
assign(objname, saom.meta)

# print results to file if requested
if (print.results) {
  saom.meta.print <- round(saom.meta, 4)
  write.xlsx(saom.meta.print, file=paste0(objname, '.xlsx'))
}

rm(print.results, saom.meta.print, saom.meta, objname, saom.model)
rm(list=ls(pattern='saom.result'))