Preliminary
I would like to rewrite the equation that you have in a slightly different manner. Since a density matrix can be written as a matrix, we can also write it down as a linear combination of elements from a basis for the space of density matrices. We can use essentially any basis to do this, but some are preferred: most notably, the Pauli basis. For a $2$-qubit system, we use the $2$-qubit Pauli group $$^{1} \mathcal{P}^{2} = \big(\frac{1}{\sqrt{2}}\{I,X,Y,Z\}\big)^{\otimes 2};$$ the two-fold tensor product of the four Pauli matrices.
Then, we can simply expand any density matrix $\rho$ as:
\begin{equation}
\rho = \sum_{P_{i} \in \mathcal{P}^{2}}p_{i}P_{i},
\end{equation}
with $p_{i} = \mathrm{tr}\big[P_{i}\rho\big]$ the 'coefficent' of $\rho$ along $P_{i}$.
The goal of quantum state tomography is then to determine the coefficients $\{p_{i}\}$, thereby characterizing $\rho$ fully. Note that we can also try to expand $\rho$ in any other basis for the space of density matrices, we then get different coefficient $\{p_{i}\}$.
Determining the coefficients $\{p_{i}\}$
The coefficients $\{p_{i}\}$ are readily determined by the fact that a measurement of an observable $A$ on the state $\rho$ has an expectation value $\langle A,\rho \rangle = \mathrm{tr}\big[A\rho\big]$. Therefore, multiple measurements of the same observable will give information to calculate these expectation values.
The first question
The Bell states are four states that together span the $2$-qubit Hilbert space. However, the $2$-qubit density matrix space is spanned by $4^{2}$ elements; so to be able to perform the expansion of $\rho$ we need $16$ linearly independent observables. You can think of this mathematically like this:
The four Bell states $\{|\Psi_{+}\rangle,|\Psi_{-}\rangle, |\Phi_{+}\rangle, |\Phi_{-}\rangle\}$ allow us to determine the expectation values for the four observables $|\Psi_{+}\rangle\langle\Psi_{+}|,|\Psi_{-}\rangle\langle\Psi_{-}|,|\Phi_{+}\rangle\langle\Phi_{+}|$ & $|\Phi_{-}\rangle\langle\Phi_{-}|$. However, you are 'missing' the $12$ cross-products that together with these four observables span the entire density matrix space. Essentially, you are only determining the diagonal elements of the density matrix expressed in the Bell basis. To determine the off-diagonal elements you need to estimate the eigenvalues of the other observables somehow, which requires some less-than straightforward operations.
The second question
To determine the expectation values for an observable $I \otimes P$ with $P \not = I$, one can indeed 'just' measure the second qubit in the corresponding eigenbasis to $P$, and leave the other qubit unaffected (essentially tracing out this qubit). Note that if $P$ is also $I$, the expectation value (of $I \otimes I$) will be equal to $1$ by definition.
However, there is a method that can be used to infer the expectation values of these observables from other measurement outcomes. As an example, note that the observable $ZZ = Z \otimes Z$ has four eigenstates, two of which are a $+1$ eigenstate ($|00\rangle \& |11\rangle$) and two of which are a $-1$ eigenstate ($|01\rangle \& |01\rangle$). This gives the spectral decomposition for $Z \otimes Z$:
$$
Z \otimes Z = (+1)|00\rangle\langle 00| + (+1)|11\rangle\langle 11| + (-1)|01\rangle\langle 01| + (-1)|10\rangle\langle 10|
$$
Measuring the $ZZ$ observable $N_{\text{tot}}$ times will give counts$^{2}$ $N_{00}, N_{01}, N_{10}$ and $N_{11}$ which together should add up to $N_{tot}$. The expectation value for observable $ZZ$ can now be estimated:
\begin{equation}
\begin{split}
\mathrm{tr}\big[(ZZ)\rho\big] =& (+1)\mathrm{tr}\big[|00\rangle\langle00|\rho\big] + (+1)\mathrm{tr}\big[|11\rangle\langle11|\rho\big] + \\
& (-1)\mathrm{tr}\big[|01\rangle\langle01|\rho\big] + (-1)\mathrm{tr}\big[|10\rangle\langle10|\rho\big] \\
\simeq& +\frac{N_{00}}{N_{\text{tot}}} +\frac{N_{11}}{N_{\text{tot}}} - \frac{N_{01}}{N_{\text{tot}}} - \frac{N_{10}}{N_{\text{tot}}}.
\end{split}
\end{equation}
Now, the observable $I\otimes Z = IZ$ can too be expanded into its eigenspaces. The thing is, the $+1$ and $-1$ eigenstates of $Z$ (i.e. $|0\rangle$ and $|1\rangle$) are both $+1$ eigenstates of $I$. Therefore, the previously 'measured' eigenstates for $ZZ$ can be used to determine the expectation value for $IZ$ as well. The states $|00\rangle$ and $|10\rangle$ are now the $+1$ eigenstates, whereas $|01\rangle$ and $|11\rangle$ are now the $-1$ eigenstates; we can thus approximate the expectation value for $IZ$ as:
\begin{equation}
\mathrm{tr}\big[IZ\rho\big] \simeq +\frac{N_{00}}{N_{\text{tot}}} +\frac{N_{10}}{N_{\text{tot}}} - \frac{N_{01}}{N_{\text{tot}}} - \frac{N_{11}}{N_{\text{tot}}}.
\end{equation}
The expectation value for $ZI$ can be determined in a similar manner from the same measurement outcomes, keeping track of the new eigenvalues for all states.
Since all states are $+1$ eigenstates of $I$, and therefore also those of the $X$ and $Y$ operators, we can use a similar technique to determine $IX, XI, IY \ \& \ YI$. This means that instead of performing $4^{2}$ different measurements $N_{\mathrm{tot}}$ times, we only need to perform $3^{n}$ different measurements. This technique also works for a higher number of qubits, which brings down the cost of QST somewhat.
Final notes
$^{1}$ Please note that I have omitted the term $\frac{1}{2^{n}}$ from the expansion of $\rho$. This means that I have implicitly assumed the Pauli basis to be the normalized Pauli basis: $\mathrm{tr}\big[P_{i}P_{j}\big] = \delta_{ij}$ (with $\delta_{ij}$ the Kronecker delta), instead of $2^{n}\delta_{ij}$. If you wish to use the normal Pauli basis, you have to include the factor of $\frac{1}{2^{n}}$ in the expansion.
Note that using a normalized basis does help in some occasions (for instance, here it clears up the notation a bit.) However, you also lose some important properties, most notably that these normalized matrices do not form a group under multiplication anymore.
$^{2}$ This is a subtle point. How is the measurement according to this observable actually implemented? If you have some magical measurement device that actually only tells you if the measured state was (projected to) a $+1$ or $-1$ eigenstate of $Z\otimes Z$, you do not learn those individual counts ($N_{00}$ etc.); but these results would still be enough to approximate $\mathrm{tr}\big[(Z\otimes Z)\rho\big]$ though! However, in practice, this measurement is emulated by measuring $Z$ on both qubits separately, which would give you counts for all four outcomes.
