In neuroimaging studies, multivariate methods provide a framework for studying associations between complex patterns distributed throughout the brain and neurological, psychiatric, and behavioral phenotypes. However, mitigating the influence of nuisance variables, such as confounders, remains a critical challenge in multivariate pattern analysis (MVPA). In studies of Alzheimer's Disease, for example, imbalance in disease rates across age and sex may make it difficult to distinguish between structural patterns in the brain (as measured by neuroimaging scans) attributable to disease progression and those characteristic of typical human aging or sex differences. Concerningly, when not properly adjusted for, nuisance variables can obscure interpretations and preclude the generalizability of findings from neuroimaging studies. Motivated by this critical issue, in this work we examine the impact of nuisance variables on features extracted from image decomposition methods and propose Penalized Decomposition Using Residuals (PeDecURe), a new MVPA method for obtaining nuisance variable-adjusted features. PeDecURe estimates primary directions of variation which maximize covariance between residualized imaging features and a variable of interest (e.g., Alzheimer's diagnosis) while simultaneously mitigating the influence of nuisance variation through a penalty on the covariance between residualized imaging features and those variables. Using features estimated using PeDecURe's first direction of variation, we train an accurate and generalizable predictive model, as evidenced by its robustness in testing samples with different underlying nuisance variable distributions. We compare PeDecURe to commonly used decomposition methods (principal component analysis (PCA) and partial least squares) as well as a confounder-adjusted variation of PCA. We find that features derived from PeDecURe offer greater accuracy and generalizability and lower partial correlations with nuisance variables compared with the other methods. While PeDecURe is primarily motivated by MVPA in the context of neuroimaging, it is broadly applicable to datasets where the dimensionality or complexity of the covariance structure calls for novel methods to handle sources of nuisance variation.