Using the latest stellar evolution models, theoretical stellar spectra, and a compilation of observed emission line strengths from Wolf-Rayet (WR) stars, we construct evolutionary synthesis models for young starbursts. We explicitly distinguish between the various WR subtypes (WN, WC, WO), whose relative frequency is a strong function of metallicity, and we treat O and Of stars separately. We calculate the numbers of O and WR stars produced during a starburst and provide detailed predictions of UV and optical emission line strengths for both the WR stellar lines and the major nebular hydrogen and helium emission lines, as a function of several input parameters related to the starburst episode. We also derive the theoretical frequency of WR-rich starbursts. Our models predict that nebular HeII lambda 4686 emission from a low-metallicity starburst should be associated with the presence of WC/WO stars and/or hot WN stars evolving to become WC/WO stars. In addition, WR stars contribute to broad components beneath the nebular Balmer lines; the broad WR component may constitute several percent of the total flux in the line. We review the various techniques used to derive the WR and O star content from integrated spectra, assess their accuracy, and propose two new formulae to estimate the WR/O number ratio from UV or optical spectra. We also explore the implications of the formation of WR stars through mass transfer in close binary systems in instantaneous bursts. While the formation of WR stars through Roche lobe overflow prolongs the WR dominated phase, there are clear observational signatures which allow the phases in which WR stars are formed predominantly through the single or the binary star channels to be distinguished. In particular at low metallicities, when massive close binaries contribute significantly to the formation of WR stars, the binary-dominated phase is expected to occur at ages corresponding to relatively low H_beta equivalent widths. The observational features predicted by our models allow a detailed quantitative determination of the massive star population in a starburst region (particularly in so-called ``WR galaxies'') from its integrated spectrum and provide a means of deriving the burst properties (e.g., duration, age) and the parameters of the initial mass function of young starbursts. The model predictions should provide the most reliable determinations to date. They can also be used to test the current theories of massive star evolution and atmospheres and investigate the variation in stellar properties with metallicity.