This study presents a mixed finite element formulation for the stress analysis of laminated composite beams based on the refined zigzag theory. The Hellinger-Reissner variational principle is employed to obtain the first variation of the functional that is expressed in terms of displacements and stress resultants. Due to C0 continuity requirements of the formulation, linear shape functions are adopted to discretize the straight beam domain with two-noded finite elements. The proposed formulation is shear locking free from nature since it introduces displacement and stress resultant terms as independent field variables. A monolithic solution of the global finite element equations is preferred, hence the stress resultants are directly obtained from the solution of these equations. The in-plane strain measures of the beam are obtained directly at the nodes over the compliance matrix and stress resultants by avoiding error-prone spatial derivatives. Following, transverse shear stresses are calculated from the equilibrium equations at the post-processing level. This simple but effective finite element formulation is first verified and tested for convergence behavior. The robustness of the approach is shown through some examples and its accuracy in predicting the displacement and stress components is revealed.