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.
More mathemtically
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.
There is no answer to this. When you are taught to use dimensional analysis at school the teacher invariably selects an easy example (it's almost always the pendulum) to keep things simple. In the real world there is no guarantee that you have a dimensionless constant.
It's actually quite rare to use dimensional analysis to derive equations in the real world. The sorts of simple systems that are amenable to dimensional analysis are usually already well known. However it's very, very useful to use dimensional analysis to check that an equation you derive is dimensionally consistent.
For example suppose you're working through a differential equation for some quantity, and after covering many sheets of paper with scribbles you end up with a final equation. It's very easy to make a minor mistake along the way, so the first thing you check is that your final equation is dimensionally consistent, i.e. the dimensions of the left and right sides are the same. If they aren't that means you've made a mistake somewhere. I routinely do this in my answers to questions on this site!
Best Answer
A major limitation of dimensional analysis is that you must know some of the physics behind the concept that you are attempting to analyze, and since you have chosen to make the proportionality constant dimensionless, it has skewed your results.
However, we can easily deduce the dimensional properties of this constant by rearranging the equation you were attempting to analyze :
$F= \dfrac{K * m_1*m_2}{r^{2}}$
$K= \dfrac{r^2*F}{ m_1*m_2}=\dfrac{m*v*r^2}{ m_1*m_2*t}$
In dimensional form:$\implies \dfrac{[M][L^1][T^{-1}L^{-2}]}{[T][M][M]}={[M^{-1}][L^3][T}^{-2}]$, which is exactly what we would expect for the dimensions of the gravitational proportionality constant. However, to get to this point, you still have to know that $F=m*a$.