Magnetorheological elastomers (MREs) are composite materials whose mechanical properties can rapidly and reversibly change with the application of a magnetic field. MREs are composed of a non-magnetic polymer matrix embedded with magnetic particles. These materials can also demonstrate magnetostriction, with large strains compared to the magnetostriction seen in crystalline solids. Ultrasoft MREs, with Young's moduli on the order of kPa or less, have characteristic pinched magnetic hysteresis loop shapes that are associated with rearrangement and clustering of magnetic particles within the surrounding polymer matrix. MREs have potential applications as vibration isolators and absorbers, soft actuators, sensors, and applications in biomedical fields. Modeling of these systems at different length and time scales is an important challenge for understanding their complex magnetomechanical behavior, fabricating MREs with desired properties, and designing MRE devices. This dissertation covers efforts centered on creating models of MRE systems at length scales that resolve the micron sized magnetic particles, but not the polymer chains of the surrounding matrix, in order to study the magnetic response, particle motion, polymer deformation, and field-dependent stiffness. A simple dipole-spring model is developed to investigate the interplay between magnetic and elastic forces in a two particle system. Results from this model show particle motion on the scale of the particle size, the role of particle clustering in magnetic hysteresis in these systems, and a pinched magnetic hysteresis loop shape that qualitatively matches experimentally measured hysteresis loops for ultrasoft MREs. Magnetometry measurements of MRE samples with polymer matrices with varying Young's moduli shows decreasing magnetic hysteresis with increasing polymer stiffness. Similarly, cooling an ultrasoft MRE below a glassy phase transition increases the elastic modulus, and magnetometry measurements show no magnetic hysteresis when below this temperature. A more sophisticated model is developed that takes into account the nonlinear magnetic susceptibility and gradual approach to saturation of the carbonyl iron particles, and treats the surrounding polymer matrix as a linearly elastic continuum that surrounds the particles. Results for two particle systems in this model demonstrate non-zero remanent magnetization and switching fields for the magnetization direction. This more sophisticated modeling approach was used to study the combined response of MREs to magnetic fields and strains. Shear strain simulations for two particle systems show the effective shear modulus increases with applied magnetic field up to magnetic saturation, that the relative size of this effect increases with the volume fraction of magnetic particles up to a point, and that particle clustering enhances the stiffening effect. Shear strain simulations for chain-like arrangements of magnetic particles were also done to study the effect of magnetic field orientation on the effective shear modulus. These simulations show that the shear modulus increases with increasing magnetic field when the field is applied parallel to the chain, but the shear modulus decreases with increasing field when the field is applied perpendicular to the chain. Finally, helical arrangements of particles with varying helix radii show the complex relationship between particle geometry and the effective shear modulus, as well as the ability for applied strains to create or break particle clusters.