Magnetorheological elastomers (MREs) are composite materials consisting of magnetizable particles embedded in an elastomeric matrix material. They are capable of magnetostriction, generating actuation traction, and magnetic field-dependent modulus effects. Because of these properties, MREs have a myriad of potential applications including magnetic position sensors, electromagnetic shielding, and flexible magnets, as well as controllable mounts, clutches and vibration absorbers. While experimental results demonstrate the promise of these materials, the effects that can be obtained are still relatively small. The goal of this thesis is to provide a better understanding of the properties of MREs using theoretical methods to help guide their continued development. For this purpose, we use homogenization, which determines an effective macroscopic constitutive model for an MRE based on the properties of the constituent phases and their arrangement within the composite. Variational homogenization methods were developed in this work which provide a framework to predict the behavior of general magnetoelastic composites. However, specializing this work to MREs, a somewhat simplified approach is developed which assumes that the microstructure evolves exactly as it would in the purely mechanical problem; we refer to it as the &ldquopartial decoupling approximation.&rdquo Specific constitutive models for MREs made with rigid inclusions are derived which incorporate the non-linear effects of magnetic saturation and the non-linearity inherent in finite strain mechanics. While the magnetoelastic coupling in MREs can be accounted for by considering the torques and forces exerted on particles by the applied magnetic field, the variational approach used here circumvents the need to explicitly compute these forces and torques. The results demonstrate that for aligned loading, where the magnetic torques vanish, the magnetoelastic coupling is proportional to the square of the particle concentrations to leading order. For non-aligned loading, the associated torques have effects proportional to the concentration and can be significantly larger. We optimize magnetoelastic properties like magnetostriction, actuation traction, and magnetoelastic modulus with respect to the microstructure. Furthermore, we investigate multi-scale composites that utilize magnetic torques and particle rotations to produce strong magnetoelastic coupling.