Why is $langle textbf R rangle =Tr rho R$ here $rho$ is the state density matrix

I believe the language in the text you quote is quite confusing and unusual (but this may only be my perception due to the books I use).

Consider an ensemble of systems, it is prepared in a way, that a system in the ensemble is in state $left| n right>$ with probability $p_n$.

A density matrix is an operator of the form: $$ rho = sum_n left|nright> p_n left<nright|. $$

Where $left| n right>$ is a complete (not necessarily orthogonal) set of states.

If we now consider an observable $R$, and consider the expression (for simplicity let $left| n right>$ be orthogonal, but it also works out otherwise): $$ text{Tr}(rho R) = sum_{m, n} left<mmiddle|n right> p_n left<nmiddle| R middle| m right> = sum_n p_n left<n middle| R middle|nright> = sum_n p_n left< R right>_n = left< R right>$$ Where we evaluate the trace in the same basis set as chosen to specify the statistical operator. That is, the considered expression evaluates to what we expect of an average (the quantum mechanical average of $R$ in the states averaged over the probability of the states)!

Observables are given by self adjoint operators, what we seek is the expectation values for these observables. So we need a probability measure that assigns to each eigenvalue it's probability of occurrence. This reduces the problem to finding the measures on separable Hilbert space $mathcal{H}$.

Spectral theorem: If $A$ be a self adjoint operator then it can be decomposed as, $$A=sum_i lambda_i P_i$$ where $lambda_i$s are eigenvalues and $P_i$ are the corresponding projection operators.

Gleason's theorem: The only possible measures on complex separable Hilbert spaces $mathcal{H}$ of dimension greater than 2 are measures of the form $$mu(lambda_i)=Tr(rho E_i).$$ where $rho$ is a positive semidefinite self adjoint operator of unit trace and $E$ is the spectral measure of the self adjoint operator $A$.

Now the expectation value of the observable will be given by, $$langle Arangle= sum_ilambda_i mu(lambda_i)=Tr(rho A)$$