  We show that short-term perturbations among massive planets in multiple
planet systems can result in radial velocity variations of the central star
which differ substantially from velocity variations derived assuming the
planets are executing independent Keplerian motions. We discuss two alternate
fitting methods which can lead to an improved dynamical description of multiple
planet systems. In the first method, the osculating orbital elements are
determined via a Levenberg-Marquardt minimization scheme driving an N-body
integrator. The second method is an improved analytic model in which orbital
elements are allowed to vary according to a simple model for resonant
interactions between the planets. Both of these methods can determine the true
masses for the planets by eliminating the sin(i) degeneracy inherent in fits
that assume independent Keplerian motions. We apply our fitting methods to the
GJ876 radial velocity data (Marcy et al. 2001), and argue that the mass factors
for the two planets are likely in the 1.25-2.0 range
