  For $m$ number of bosons, carrying spin ($\cs=\spin$) degree of freedom, in
$\Omega$ number of single particle orbitals, each doubly degenerate, we
introduce and analyze embedded Gaussian orthogonal ensemble of random matrices
generated by random two-body interactions that are spin ($S$) scalar
[BEGOE(2)-$\cs$]. Embedding algebra for the BEGOE(2)-$\cs$ ensemble and also
for BEGOE(1+2)-$\cs$ that includes the mean-field one-body part is $U(2\Omega)
\supset U(\Omega) \otimes SU(2)$ with SU(2) generating spin. A method for
constructing the ensembles in fixed-($m,S$) spaces has been developed.
Numerical calculations show that for BEGOE(2)-$\cs$, the fixed-$(m,S)$ density
of states is close to Gaussian and level fluctuations follow GOE in the dense
limit. For BEGOE(1+2)-$\cs$, generically there is Poisson to GOE transition in
level fluctuations as the interaction strength (measured in the units of the
average spacing of the single particle levels defining the mean-field) is
increased. The interaction strength needed for the onset of the transition is
found to decrease with increasing $S$. Covariances in energy centroids and
spectral variances are analyzed. Propagation formula is derived for the
variance propagator for the fixed-$(m,S)$ ensemble averaged spectral variances.
Variance propagator clearly shows, by applying the Jacquod and Stone
prescription, that the BEGOE(2)-$\cs$ ensemble generates ground states with
spin $S=S_{max}$. This is further corroborated by analyzing the structure of
the ground states in the presence of the exchange interaction $\hat{S}^2$ in
BEGOE(1+2)-$\cs$. Natural spin ordering is also observed with random
interactions. Going beyond these, we also introduce pairing symmetry in the
space defined by BEGOE(2)-$\cs$. Expectation values of the pairing Hamiltonian
show that random interactions exhibit pairing correlations in the ground state
region.
