The other parts, other than the inverse square, were clear already before Newton, or at least were easy to guess. That the force of gravity is proportionality to mass of a small object responding to the field of another comes from Galileo's observation of the universal acceleration of free fall. If the acceleration is constant, the force is proportional to the mass. By Newton's third law, the force is equal and opposite on the two objects, so you can conclude that it should be proportional to the second mass too.
The model which gives you this is if you assume that everything is made from some kind of universal atom, and this atom feels an inverse square attraction of some magnitude. If you sum over all the pairwise attractions in two bodies, you get an attraction which is proportional to the number of atoms in body one times the number of atoms in body two.
So the only part that was not determined by simple considerations like this was the falloff rate. I should point out that if you look at two sources of a scalar field, and look at the force, it is always proportional to $g_1$ times $g_2$, where $g_1$ and $g_2$ are the propensity of each source to make a field by itself. Further, if you put two noninteracting sources next to each other, this g is additive, if the field is noninteracting, essentially for the reasons described above--- the independent attractions are independent. So that the proportionality to an additive body constant you multiply over the two bodies is clear. That for gravity, the g is the mass, this was established by Galileo.
Let's call the force law between the objects $F(m_1,m_2,r)$. We know that if we put the body m_1 in free fall, the acceleration doesn't depend on the mass, so
$$ F(m_1,m_2,r) = m_1 G(m_2,r) $$
So that the mass will cancel in Newton's law to give a universal acceleration. This gives you the relation
$$ F( a m_1 , m_2, r ) = a F(m_1,m_2,r) $$
We know that if we put body 2 in free fall, the same cancellation happens, but we also know Newton's third law: $F(m_1,m_2,r)= F(m_2,m_1,r)$ so that
$$ F( m_1, a m_2, r) = a F(m_1,m_2,r) $$
So you now write
$$ F( m_1 \times 1 , m_2 \times 1, r) = m_1 F( 1, m_2\times 1 , r) = m_1 m_2 F(1,1,r) $$
And this tells you that the force is proportional to the masses times a function of r. The form of the function is undetermined.
An independent argument for the scaling is that if you consider the object m_1 as composed of two nearby independent objects of mass $m_1/2$, then
$$ F(m_1/2 , m_2 , r) + F(m_1/2 , m_2 , r) = F(m_1,m_2,r)$$
Then the same conclusion follows.
These types of scaling arguments are second nature by now, and they are automatically done by matching units. So if you have a force per unit mass, the force between two massive particles must be per unit mass 1 and per unit mass 2.
This general argument fails for direct three-body forces, where the force between 3 bodies is not decomposable as a sum of forces between the pairs bodies individually. There are no macroscopic examples, since the pairwise additivity is true for linear fields, but the force between nucleons has a 3-body component.