Radiomic approaches in precision medicine are promising, but variation associated with image acquisition factors (“batch effects”) can affect study reproducibility and model generalizability, preventing the clinical translation of radiomic predictive models that could otherwise improve the quality of patient care. While approaches exist for the detection and correction of batch effects, they are unable to effectively detect batch effects in non-Gaussian data or simultaneously harmonize by multiple batch effects. In this work, I use PERMANOVA (Permutational Analysis of Variance) testing and RESI (Robust Effect Size Index) quantification to improve batch effect detection over traditional statistical testing in the settings of non-Gaussian data and large sample sizes, respectively. To address the problem of harmonization by multiple batch variables, I developed two approaches. The first, called OPNested ComBat, iteratively applies ComBat and identifies the optimal order of batch variables to achieve harmonization by multiple imaging parameters. The second, called MultiComBat, is a novel estimation methodology that modifies the method of moments estimation in the ComBat algorithm to enable simultaneous estimation of the parameters needed to correct for multiple imaging acquisition parameters. The results of this work will make available improved batch effect detection and harmonization methods, enabling greater study reproducibility and increasing the clinical translation of imaging-based feature applications.