What does the sum of two qubits tell about their correlations? 
How much can I learn about correlations between two quits by measuring
  the sum of their values? What is the best way to formalize such a
  question?

Below is my original, longer formulation of the same question in spin language.
I have an ensemble of pairs of spin-$1/2$ particles. All I can do is to measure total $S_z$ for each pair, each time learning  one of the three standard outcomes: $+1$, $0$ or $-1$. Based on the knowledge of this probability distribution $P(S_z)$, what can be said about the density matrix of my ensemble? In particular, what is best measure for correlations between the two particles in these circumstances? 
I understand that I won't be able to distinguish whether the correlations are classical or necessarily quantum-mechanical, but struggle to come up with a measure better than a classical correlation coefficient $$r=\frac{\langle s_1 s_2 \rangle-\langle s_1 \rangle \langle s_2 \rangle}{\sqrt{\langle (s_1- \langle s_1 \rangle)^2 \rangle \langle (s_2- \langle s_2 \rangle)^2 \rangle }}$$  where $s_1$ and $s_2$ are the $\pm 1/2$ spin projections (or $0$, $1$ qubit values) of the 1st and the 2nd particle, respectively, so that my only observable can be expressed as $S_z=s_1+s_2$. I then minimize $r$ over all possible   $\langle s_1 \rangle$ and $\langle s_1 \rangle$ consistent with measured $P(S_z)$ and declare that the outcome is the minimal possible correlation in the pair. Is there a more elegant way?
UPDATE:
Classically, if we take the definition of zero correlations as $\langle s_1 s_2 \rangle =\langle s_1  \rangle \langle s_2 \rangle$, then it implies that $s_1 = +1/2$ and $s_2=+1/2$ with some independent probabilities $p_1$ and $p_2$, so that the sum is distributed bimomially $P(+1)=p_1 p_2$ and $P(-1)=(1-p_1)(1-p_2)$. These two equations have a solution with $0\leq p_1, p_2 \leq 1$ only if $$P(+1) < 1-2 \sqrt{P(-1)}+P(-1)$$
Graphically, this condition is satisfied within the region bounded by a black line in the plot below. The rest of the available probability space is the region of excess correlations ($r_{\text{min}}>0$, depicted in red). 

Updated question: is there a quantum loophole in such analysis? I suspect not but would really like to come up with a shorter/cleaner argument. 
Update no. 2 My selected spin/qubit degree of freedom is one of the many variables characterizing the state of the two particles. Can we be sure there exists no pure product state that can lead to $r_{\text{min}}>0$? See a related active question on fermionic entanglement of two particles.
Update no. 3 The work that motivated this question has been published in Nature Nanotechnology; the open access version is arXiv:1404.0030.
 A: Edited to add the second part
Edited again, for part 3 and 4
$\newcommand\ket[1]{\left|#1\right>} \newcommand\bra[1]{\left<#1\right|} $
1. Absence of Quantum Loophole
You can easily see that there is no "quantum loophole" in your argument by writing explicitly any pure separable state. With your notations, we have :
$$ 
\left(\sqrt{p_1}\ket\uparrow+e^{i\phi_1}\sqrt{1-p_1}\ket{\downarrow}\right)
\left(\sqrt{p_2}\ket\uparrow+e^{i\phi_2}\sqrt{1-p_2}\ket{\downarrow}\right)
$$
This expression obviously corresponds to all possible separable pure state and gives the uncorrelated distribution you describe.
Of course, as @Norbert Schuch said, if you allow for mixed states, you have access to whole region, which is the same as saying that you allow classical correlations.
To detect quantum correlations, as you are probably aware, you will need to have information about another direction, like $S_x$.
2. Lower bound of entanglement of a pure state
If you are sure that your state $\ket\psi$ is pure, you can have a lower bound of its entanglement. Let $P_\pm$ and $P_0$ be the three spin probabilities. The state you have can be written 
$$\ket{\psi}=\sqrt{P_+}\ket{\uparrow\uparrow} + e^{i\phi}\sqrt{P_-}\ket{\downarrow\downarrow} +\sqrt{P_0}(\alpha\ket{\uparrow\downarrow}+\beta\ket{\downarrow\uparrow}) \text{ with } |\alpha|^2+|\beta|^2=1$$
Since our measurement does not distinguish what happens in the $S_Z=0$ subspace, we probably have no information about the entanglement there. So we can take the fidelity $F$ with a maximally entangled correlated state for entanglement measure. We have
$$
F:=\max_{\theta}\left|\left<\psi\right|\frac{\ket{\uparrow\uparrow}+e^{i\theta}\ket{\downarrow\downarrow}}{\sqrt2}\right|^2=\frac{\left(\sqrt{P_+}+\sqrt{P_-}\right)^2}2
$$
You cannot prove entanglement when $F=\tfrac12$. This corresponds to the zero correlation line :-)
Of course, all this is only valid for pure states, and could also correspond to classically correlated states.
3. Lower bound, take 2
Since the state is pure, its entanglement is fully characterized by the eigenvalues $\lambda_\pm$
 of its partial trace $\rho=Tr_B \ket\psi \bra\psi$. They are the solutions of 
$\lambda^2-\lambda+\det\rho$ and are given by
$\lambda_\pm=\tfrac12 \pm \sqrt{\tfrac14 -\det\rho}.$
We have the Schmidt decomposition
$\ket{\psi}=\sqrt{\lambda_+}\ket{0}\otimes\ket{0}+\sqrt{\lambda_-}\ket{1}\otimes\ket{1}$,
where the $\ket{i}$s are the bases which diagonalise $ Tr_B \ket\psi \bra\psi$ and 
$ Tr_A \ket\psi \bra\psi$.
Before computing $\det\rho$, one can see that $0\le\det\rho\le\tfrac14$. When $\det\rho=0$, $\lambda_+= 1$, $\lambda_-=0$ and
$\ket\psi$ is separable ; when $\det\rho=\tfrac14$, 
$\lambda_+=\lambda_-=\tfrac12$ and
$\ket\psi$ is maximally entangled. 
The bigger $\det\rho$ is, the more $\ket\psi$ is entangled.
More quantitatively, the entropy of entanglement of this state is $h(\lambda_-)$ where $h(\cdot)$ is the binary entropy function and the negativity of the state is $\det\rho$.
\begin{align}
\det\rho
  &=\left(P_++P_0|\alpha|^2\right)\left(P_-+P_0|\beta|^2\right)
  -\left| \beta\sqrt{P_0P_+}+\alpha^*e^{i\phi}\sqrt{P_0P_-}\right|^2
  \\
 &=P_+P_- +|P_0\alpha\beta|^2 +P_0\left(
   |\alpha|^2P_-+ |\beta|^2P_+
   - |\beta|^2P_ +- |\alpha|^2P_-
  - 2 Re\left( \beta \alpha^*e^{i\phi}\sqrt{P_+P_-}
  \right)\right)\\
  &=\left| \sqrt{P_+P_-}   + \beta \alpha^*e^{i\phi} P_0 \right|^2
\end{align}
We want a lower bound on the entanglement. We therefore chose the free parameters
$\alpha$, $\beta$, and $\phi$ in order to minimize $\det\rho$. Since we have 
$|\alpha\beta|\le\tfrac12$, $\min\det\rho=0$ if $ \sqrt{P_+P_-}\le \tfrac12 P_0$. 
This corresponds to the separability limit established earlier.
If we have 
$|\alpha\beta|\le\tfrac12$, then
$$
  \min\det\rho=\left( \sqrt{P_+P_-}- \tfrac12 P_0\right)^2
   =\tfrac14 \left[\left(\sqrt{ P_+}+ \sqrt{ P_-}\right)^2 - 1\right]^2
   = \left(F-\tfrac12\right)^2.
$$
For these minimally entangled states, $\lambda_\pm=\frac12\pm\sqrt{F(1-F)}$.
4. Bonus : an upper bound on the entanglement
One can also look at an upper bound on the entanglement with the same approach:
$$
  \max\det\rho= \left| \sqrt{P_+P_-}   + \tfrac12 P_0 \right|^2
    = \left(G-\tfrac12\right)^2 \text{ with } G=\tfrac12\left( \sqrt{ P_+} - \sqrt{ P_-}\right)^2.
$$
For these "maximally" entangled states, $\lambda_\pm=\frac12\pm\sqrt{G(1-G)}$.
A good way to put all this on a graph is to take $\sqrt{P_\pm}$ as coordinates. The possible measurement are in the top right quarter of a disk and the bound described here are constant along 45° straight lines. 
A: Computation of the $S_z$ probability distribution for each of the
manifolds of equal entanglement: 
Remark: Notations and references from  Kuś and Žyczkowski are used.
Case 1: The separable case: The state vector is parametrized as (equation: 24)
$w = \begin{bmatrix} \cos \alpha \cos \beta e^{i \chi_1},& \cos \alpha \sin\beta  e^{i \chi_2}, &\sin \alpha \cos\beta e^{-i \chi_1}, & \sin \alpha \sin \beta e^{-i \chi_1} \end{bmatrix}^{t}$
Here is the computation of the the probability distribution of $S_z$ for
the most general pure two qubit state vector:
The total spin projector operator is given by
$S_z = \frac{1}{2} \sigma_3 \otimes 1_2 +   1_2 \otimes \frac{1}{2}\sigma_3$
The characteristic function of the probability distribution is given by:
$ \Phi_{S_z}(t) = w^{\dagger} e^{i t S_z} w  =P(S_z = 1) e^{it} + P(S_z = 0)+ P(S_z = -1) e^{-it}$
From which the probability distribution can be read:
$P(S_z = 1) =  \cos^2 \alpha\cos^2 \beta $
$P(S_z = 0) = \cos^2 \alpha\sin^2 \beta+ \sin^2 \alpha\cos^2 \beta $
$P(S_z = -1) =  \sin^2 \alpha\sin^2 \beta $
Case 2: The maximally entangled case: The state vector is parameterized as (equation: 23)
$w = \frac{1}{\sqrt{2}}\begin{bmatrix} \cos \alpha  e^{i \chi_1},& \sin \alpha   e^{i \chi_2}, &\sin \alpha  e^{-i \chi_1}, & \cos \alpha  e^{-i \chi_1} \end{bmatrix}^{t}$
In this case:
$P(S_z = 1) =  \frac{1}{2}\cos^2 \alpha$
$P(S_z = 0) =  \sin^2 \alpha$
$P(S_z = -1) =  \frac{1}{2}\cos^2 \alpha $
Case 3: The generic case: The state vector has been worked out:
$w = \frac{1}{\sqrt{1+\sin^2\eta}} \begin{bmatrix} (\cos \alpha \cos \beta -\sin\eta e^{i\phi}\sin\alpha \sin\beta) e^{i \chi_1},& (\cos \alpha \sin\beta + \sin\eta e^{i\phi} \sin\alpha \cos \beta) e^{i \chi_2}, \\&(\sin \alpha \cos\beta + \sin\eta e^{i\phi} \cos\alpha \sin \beta) e^{-i \chi_1}, & (\sin \alpha \sin \beta  -\sin\eta e^{i\phi} \cos\alpha\cos\beta) e^{-i \chi_2} \end{bmatrix}^{t}$
Please, observe that this state vector interpolates between the separable case at $\eta =0$ and the maximally entangled case at $\eta =\frac{\pi}{2}  $ (in this case it depends  on a single combination $\alpha ,  \beta, \phi$  making two of them redundant.
This state vector describes the generic orbit for $0<\eta< \frac{\pi}{2}$.
The angle $\eta$ is related to the angle $\theta$ used in the article by the relation:
$ \sin\theta = \frac{2\sin\eta}{\sqrt{1+\sin^2\eta}}$
In this case, we obtain:
$P(S_z = 1) = \frac{1}{1+\sin^2\eta} ( \cos^2 \alpha\cos^2 \beta + \sin^2 \alpha\sin^2 \beta \sin^2 \eta  - 2\cos\alpha\cos \beta \sin \alpha\sin \beta \sin \eta \cos \phi)$
$P(S_z = 0) = \frac{1}{1+\sin^2\eta} ( \cos^2 \alpha\sin^2 \beta (1+\sin^2\eta)+ \sin^2 \alpha\cos^2 \beta (1+\sin^2\eta)  + 4\cos\alpha\cos \beta \sin \alpha\sin \beta \sin \eta \cos \phi)$
$P(S_z = -1) = \frac{1}{1+\sin^2\eta} (\sin^2 \alpha\sin^2 \beta + \cos^2 \alpha\cos^2 \beta \sin^2 \eta  - 2\cos\alpha\cos \beta \sin \alpha\sin \beta \sin \eta \cos \phi)$
